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Preface 


The  purpose  of  this  study  was  to  conduct  an  analytical 
experiment  in  the  application  of  Floquet  theory  to  the  heli¬ 
copter  rotor  blade  flapping  motion  problem.  Prior  to  this 
time,  Floquet  application  had  been  restricted  to  examining 
the  stability  of  the  blade  flapping  motion.  This  study  extends 
the  application  to  using  Floquet  theory  as  a  tool  for  blade 
feedback  control  design  synthesis. 

In  conducting  this  study,  several  computer  algorithms 
were  used  to  derive  the  key  information  from  the  blade  equations 
of  motion  needed  to  design  the  feedback  controller.  All  the 
computer  software  needed  to  conduct  this  study  is  included  in 
the  appendices  of  the  text  in  their  order  of  implimentation. 

This  software  is  written  in  Texas  Instruments  Extended  Basic, 
but  is  easily  convertible  to  another  computer  language.  Out¬ 
put  from  the  algorithms  is  also  included  in  the  appendices 
for  model  comparisons . 

In  conducting  this  analysis,  I  am  deeply  indebted  to 
others  for  their  encouragement  and  assistance.  First,  a  word 
of  thanks  to  my  thesis  advisor,  Dr.  Robert  Calico,  and  to  Dr. 
William  Wiesel  for  the  countless  hours  of  discussions  and  in¬ 
sights  received  on  the  idiosyncrasies  of  Floquet  theory  which 
allowed  me  to  apply  it  to  this  study.  And,  a  sincere  thanks 
to  my  wife  Jane  and  daughters  Elizabeth  and  Victoria  for 
their  understanding  and  unending  support  through  this  some¬ 
times  trying  odyssey. 


James  K.  March 


Table  of  Contents 


Preface  . . 

List  of  figures  . 

List  of  Tables  . . . 

List  of  Symbols  . 

Abstract  . 

I .  Introduction  . 

II.  Helicopter  Rotor  Flapping  Dynamics  .... 

Coordinate  Frame  . 

Flapping  Hinge  Moments  . 

Blade  Aerodynamics  . . 

Stability  Model  . 

III .  Floquet  Theory  . 

Foundations  of  Floquet  Theory  . . . . 
Feedback  Control  Synthesis  . 

Single  Mode  . 

Single  Scalar  Modal  Control 
Scalar  Control  of  Two  or  More 
Modes  . 

Recapitulation  . 

IV.  Results  . . 

No  Feedback  Case  . 

Feedback  Case  . .  .  . . 

V.  Conclusion  . 

Appendix  As  Numerical  Methods  . 

Integration  Routines  . 

Fourier  Series  Coefficients 

Analysis  Routines  . 

Eigenproblem  Routine  . 

Feedback  Routine  . . 

Step  Size  Determination 
Analysis  Procedure  . 


Appendix  B:  F  Matrix  for  Y  =8  66 

Appendix  C:  Fourier  Series  Coefficients 

of  F  Matrices  .  73 

For  X  =2  74 

For  r  =4  78 

For  y  =6  82 

For  y  =8  86 

For  y  =10  90 

For  y  =12  94 

For  y  =14  98 

For  y  =16  102 

Appendix  D:  Mode  Controllability  Matrices 

Element  Fourier  Series 

Coefficients  .  106 

For  y  =2  10? 

For  y=4  Ill 

For  y  =6  115 

For  y  =8  119 

For  y  =10  123 

For  y  =12  127 

For  y  =14  131 

For  y=l6  .  135 

Appendix  E«  Monodromy  Matrix  Routine  .  139 

Appendix  F:  Eigenvalue/Eigenvector  Routine  .  144 

Appendix  G:  F  Matrix  Routine  . 165 

Appendix  H:  Fourier  Coefficients  Analysis 

Routine  .  173 

Appendix  I:  Mode  Controllability  Matrix 

Routine  .  176 

Appendix  J:  Feedback  Controller  Routine  .  I83 

Bibliography  .  194 

Vita  .  198 


Figure 

1 .  Blade  Coordinate  Frame  . . . .  . . 

2.  Blade  External  Forces  . . 

3 .  Blade  Aerodynamics  . 

4.  Flow  Regions  . . . . 

5.  Stable  and  Unstable  Response  ... 

6.  Stability  Boundary  . . . 

7.  Root  Locus  . 

8.  Feedback  Control  Performance  ... 


List  of  Symbols 


lift  coefficient  versus  angle  of 
attack  curve  slope 

Fourier  series  coefficients 


tip  loss  factor 

aerodynamic  damper  characteristic 
of  blade  on  flap  moment 

chord  of  airfoil  section 

lift  coefficient  for  a  given  blade 
element  angle  of  attack 

differential  centrifugal  force 

forcing  function  due  to  aerodynamic 
forces  on  blade 

Fourier  series 

gravitational  acceleration  constant 

mode  controllability  matrix 

blade  moment  of  inertia  about  the 
Zhinge  axis 

aerodynamic  spring  characteristic 
on  flap  moment 

gain  matrix 

flap  restraint  spring  constant 

lift  generated  by  blade 

element  mass  of  blade  element 

blade  natural  frequency  in  hover 

total  length  of  blade 

distance  along  blade  centerline 
measured  from  hinge 

period  of  the  periodic  system 


\  *.  \ 


up 

velocity  component  in  Z.  . 
direction  ninge 

UR 

radial  velocity  component  in 

Xblade  direction 

velocity  component  in  Z.  . 
direction  ninge 

u 

forward  velocity  vector  of  heli¬ 
copter 

X 

nondimensional  distance  ratio; 
r/R 

X  Y  Z 

hinge'  hinge’  hinge 

coordinate  reference  frame  of 
hinge 

Xblade’  Yhlade ’  Zblade 

coordinate  reference  frame  of 
blade 

a 

blade  element  angle  of  attack 

0 

blade  flap  angle 

•  — 

0,0 

ls^  and  2nd  derivative  of  blade 
flap  angle  with  respect  to  time 

0,0* 

1st  and  2nd  derivative  of  blade 
flap  angle  with  respect  to  blade 
azimuth  angle 

r 

blade  lock  number 

€ 

angular  width  of  mixed  flow  region 
across  the  rotor  plane 

V 

modal  variable 

9 

blade  element  pitch  angle  to 
horizontal  Xhinge,  Yhinge  plane 

9o 

pitch  angle  of  blade  root 

9i 

built  in  linear  twist  angle  along 
blade's  longitudinal  axis 

cyclic  pitch  control  terms 

X 

inflow  ratio 

Xi 

characteristic  value 

Vlll 


blade  advance  ratio 


mass  density  in  free  stream 


state  transition  matrix 


velocity  vector  angle  to  X 
Yhinge  Plane 


hinge’ 


blade  azimuth  angle  measured  from 
tail  of  helicopter  in  counter-clock 
wise  direction 


rotor  angular  velocity 
Poincare’  exponent 


AF IT/GAE/AA/84D- 1 3 


Abstract 
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7  The  purpose  of  th  1^ xnv e s t i-gajt-ion-  was  to  explore  the 
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flapping  stability  of  a  helicopters  blade  in  forward  flight. 

The  equations  of  motion  for  the  flapping  motion  of  the  blade 
were  converted  from  nonlinear  differential  equations  with 
periodic  coefficients  to  linear  periodic  differential  equations 
through  the  assumption  of  a  rigid  blade  where  the  elastic 
flapping  deflections  are  negligible  as  compared  to  the  rigid 
body  flapping  rotations  about  the  flapping  hinge.  Aeroelastic 
effects  were  not  considered. 

The  stability  of  the  homogeneous  part  of  the  flapping 
motion  linearized  periodic  differential  equations  was  examine 
through  the  application  of  Floquet  theory.  The  flapping  blade 
motion  was  simulated  over  one  period  to  derive  the  elements 
of  the  monodromy  matrix.  The  monodromy  matrix  was  next  trans¬ 
formed  into  Jordan  normal  form  through  a  similarity  trans¬ 
formation  to  obtain  its  characteristic  values  and  eigenvectors. 
The  characteristic  values  were  converted  to  their  respected 
Poincare’  exponents  and  the  periodic  eigenvectors  composition 
was  determined  and  transformed  into  Fourier  series  represen¬ 
tations  . 

A  feedback  controller  was  constructed  using  Floquet  theory 

2 

for  the  unstable  blade  flapping  motion  case,  where;  P  =1.0, 

^  =  2.4,  7=8,  and  B=0.97*  The  performance  oP  the  feedback 
controller  was  then  compared  to  the  unstable  case. 
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APPLICATION  OF  FLOQUET  THEORY  TO  HELICOPTER 
BLADE  FLAPPING  STABILITY 

I .  Introduction 

The  helicopter  rotor  blade  may  be  modelled  as  a  spring 
mass  damper  system  with  periodic  changing  spring  and  damper 
characteristics.  The  origin  of  these  characteristics  are  due 
to  the  elastic  blade  exposed  to  varying  aerodynamic  forces  as 
well  as  aeroelastic  effects.  In  forward  flight,  these  effects 
are  more  pronounced  as  the  forward  speed  of  the  helicopter 
increases  due  to  the  variable  flow  region  across  the  rotor 
plane . 

The  first  successful  helicopters  were  made  possible 
largely  by  the  development  of  the  fully  articulated  rotor  hub, 
where  blade  motion  was  accommodated  by  flapping  and  lead-lag 
hinges  at  the  blade  root  to  reduce  the  bending  stresses  at  the 
root  (33)-  This  technique  is  similar  in  principle  to  increasing 
the  thickness  of  a  propeller  blade  for  structural  consider¬ 
ations  (43j285). 

The  root  blade  itself  can  be  considered  as  an  elastic 
beam.  The  resulting  equations  of  motion  of  the  blade  are 
nonlinear  differential  equations  with  periodic  coefficients. 

By  assuming  small  blade  pitch  and  twist  along  with  accounting 
for  the  nonlinear  terms  in  the  inertial  and  aerodynamic  forces, 
these  equations  can  be  transformed  into  linear  time -periodic 
differential  equations.  With  the  further  assumption  of  low 


forward  velocity,  they  can  be  treated  as  linear  constant 
coefficient  differential  equations. 

Initially,  the  differential  equations  were  assumed  to  have 
constant  coefficients  due  to  the  low  velocity  of  earlier  heli¬ 
copter  designs.  As  the  forward  speed  increased,  the  period¬ 
icity  became  more  apparent.  The  first  attempt  to  handle  this 
problem  converted  the  blade  equations  of  motion  to  the  form 
of  the  generalized  Hill  equation  (21).  This  approach  permitted 
the  stability  determination  of  the  blade  motion  through  the 
properties  of  the  solution  of  the  Hill  equation. 

Through  the  introduction  of  the  digital  computer,  stability 
determination  was  donducted  by  simulating  the  motion  of  the 
blade  over  several  revolutions.  However,  with  the  increase 
in  complexity  of  the  mathematical  model  due  to  the  incorporation 
of  detailed  aeroelastic  effects  on  the  blade,  the  efficient 
application  of  Floquet  theory  was  discovered.  Since  that  time, 
Floquet  theory  has  been  the  standard  technique  to  parametrically 
study  the  stability  boundaries  of  helicopter  blade  motion  for 
various  designs  and  flight  conditions . 

The  purpose  of  this  research  is  to  investigate  the  flapping 
stability  of  a  helicopter  blade  in  forward  flight.  The  first 
step  in  this  investigation  was  to  determine  the  blade  flapping 
equation  of  motion.  The  assumptions  made  in  this  development 
were  that  the  blade  was  rigid  with  no  twist  and  that  the  flap 
and  pitch  angles  were  small.  The  cross  sectional  center 
of  mass  was  also  assumed  to  coincide  with  the  elastic  axis  of 


the  blade .  Thus ,  the  only  momentum  terms  that  needed  to  be 
considered  in  the  development  were  those  due  to  inertial  forces 
aerodynamic  forces,  and  flapping  restraint  contributions. 

Other  assumptions  made  in  this  development  are  addressed 
further  in  the  text  of  this  study. 

The  next  objective  of  this  research  was  to  explore  the 
applicability  of  Floquet  theory  to  this  problem.  The  most 
obvious  application  of  this  theory  was  to  investigate  the 
stability  of  the  rotor  blade  for  different  flight  and  design 
conditions.  However,  the  use  of  Floquet  theory  was  carried 
one  step  further  by  experimenting  with  its  application  to 
feedback  control  synthesis.  This  technique  has  been  applied 
to  other  control  design  problems  dealing  with  periodic  systems 
(3,4,44). 

One  comment  which  is  made  is  that  Floquet  theory  predicts 
the  properties  and  characteristics  of  the  solution,  but  it 
does  not  generate  the  solution  itself.  To  obtain  the  solution 
to  the  Floquet  problem,  numerical  methods  are  required.  The 
numerical  methods  used  in  this  study  will  also  be  presented. 

To  understand  the  impact  of  Floquet  theory  on  feedback 
controller  design,  it  is  worth  noting  the  past  basic  design 
steps  which  have  been  used.  The  initial  blade  controllers 
consisted  of  a  controllable  pitch  anti-torque  rotor  for  control 
about  the  vertical  axis,  and  a  cyclic  pitch  control  for  the 
main  rotor  blades  for  control  of  the  other  two  axes  (45«63). 

The  gains  associated  with  these  type  of  controllers  were  deter¬ 
mined  by  first  assuming  low  forward  speed.  Through  the  use  of 


classical  and  modern  control  theory,  gain  selection  was  made 
(16,18,30,3^)-  Once  the  desired  control  properties  were  obtained 
for  this  time  invarient  case,  the  forward  speed  of  the  heli¬ 
copter  was  increased  and  the  articulated  rotor  was  digitally 
simulated  for  several  conditions  to  establish  stability  bound¬ 
aries.  The  final  gains  for  these  controllers  were  thus  deter¬ 
mined  by  varying  the  gains  until  the  poles  of  the  dynamic  region 
lie  within  the  predicted  stable  region  of  response. 

Another  approach  taken  was  to  identify  the  components  of 
blade  motion,  then  build  a  controller  to  control  these  com¬ 
ponents  using  the  same  approach  as  previously  discussed  (1,20, 
35,36).  Floquet  theory  was  then  applied  to  parametrically 
study  the  stability  boundaries  of  the  system. 

The  introduction  of  Floquet  theory  will  change  this  de¬ 
sign  process  by  giving  the  designer  a  direct  access  to  pole 
location  control  for  the  periodic  system.  Rather  than  be¬ 
ginning  at  the  time  invarient  condition  of  the  blade  flapping 
equation  of  motion,  the  designer  will  be  able  to  start  from 
the  time  periodic  condition.  This  study  will  show  why  this 
can  be  done  and  will  give  examples  of  the  control  design 
process . 

Because  this  is  a  blade  motion  stability  analysis,  a 
detailed  aeroelastic  effects  study  will  not  be  necessary. 

However,  it  should  be  mentioned  that  the  techniques  discussed 
in  this  study  can  be  equally  applied  in  the  aeroelastic' 
stability  problem  (5.10,12,13). 


ii.  flglicgjter  Rotor  Flapping  Dynasties 


As  described  in  the  previous  section,  the  helicopter  rotor 
blade  can  be  modelled  as  a  spring  mass  damper  system  with  peri¬ 
odic  changing  spring  and  damper  characteristics .  These  charac¬ 
teristics  are  excited  by  the  aerodynamic  forces  due  to  inflow 
about  the  rotor  blade,  blade  twist,  cyclic  pitch  control,  and 
collective  pitch  control  (4l).  Thus,  the  blades  are  hinged 
so  as  to  prevent  the  transmission  of  the  alternating  moments  of 
the  lateral  aerodynamic  forces  on  the  blade  in  forward  flight 
onto  the  rotor  shaft  ( 21 , 28  *  150 ) .  The  derivation  of  the 
flapping  dynamic  equation  of  motion  requires  the  knowledge  of 
these  aerodynamic  forces  and  moments,  along  with  the  elastic 
and  inertial  properties  of  the  blade. 

The  helicopter  rotor  can  be  considered  as  being  composed 
of  elastic  blades.  The  exact  description  of  the  flapping  dy¬ 
namics  requires  a  detailed  consideration  of  the  elastic  flex¬ 
ibility  of  each  blade  (38)-  There  is,  however,  a  wide  array 
of  assumptions  that  can  be  made  which  lead  to  a  variety  of 
models,  starting  with  the  simple  and  computationally  efficient 
model  to  those  which  are  capable  of  simulating  unsteady  flow 
in  high  detail  (?)•  The  decision  as  to  what  assumptions  to 
make  and  what  model  to  use  depends  on  the  purpose  of  the  model 
and  the  penalty  one  is  willing  to  pay  for  more  sophistication 
(38). 

Keeping  in  mind  that  the  purpose  of  this  study  is  to 
examine  the  blade  flapping  stability  and  control  in  forward 
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flight,  it  is  justifiable  to  assume  a  rigid  blade  where  the 
elastic  flapping  deflections  are  negligible  compared  to  the 
rig^d  body  flapping  rotations  about  the  flapping  hinge  (38). 

It  has  been  shown  (8,14,33,38)  that  this  approach  can  also  be 
used  for  hingeless  rotors .  This  is  done  by  specifying  a  rigid 
blade  with  a  selected  offset  hinge  and  an  appropriate  spring 
constant  to  match  the  properties  of  the  hingeless  rotor. 

The  purpose  of  this  section  is  to  present  the  derivation 
of  the  dynamic  flapping  equation  of  motion  of  the  blade.  This 
will  be  accomplished  using  blade  element  theory  to  represent 
the  inertial,  aerodynamic,  and  restoring  moments  on  the  blade 
flapping  motion.  Through  a  coordinate  transformation,  the 
equation  will  be  converted  to  a  usable,  state  space  represent¬ 
ation  which  was  modelled  to  conduct  the  stability  and  control 
investigation  of  this  study. 

It  is  interesting  to  note  here  why  blade  element  theory 
was  employed  in  the  development  of  these  equations .  Blade 
element  theory  has  been  extensively  used  and  has  strong  foun¬ 
dations  in  propeller  theory  (43*326).  It  allows  for  the  con¬ 
venient  use  of  experimentally  derived  airfoil  data  for  predicting 
blade  performance  to  within  10  percent  accuracy,  depending  on 
the  quality  of  airfoil  data  used. 


Coordinate  Frame 

The  coordinate  reference  frame  used  in  this  derivation 
can  be  found  in  Figure  1 .  The  Xh^nge  and  Yhinge  axes  define 
the  horizontal  plane  of  the  hinge.  Zhj_nge  defined  as  being 
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perpendicular  to  the  Xhingg  and  Yhinge  axes  and  is  positive 
upwards . 

The  Xblade  axis  is  defined  along  the  centerline  of  the 
blade.  The  offset  angle  between  the  Xfe^ade  and  xh^nge  axes 
define  the  flap  angle  /3  .  The  axis  corresponds  to  the 

Yhinge  ax^s’  whi-le  the  Zblade  ax^s  is  offset  by  the  angle  /3 

from  the  Z,  . _  axis. 

hinge 

The  blade  rotates  about  the  Z..„  axis  at  the  constant 

hinge 

angular  velocity  ft  .  The  azimuth  angle  \ft  is  measured  from  the 
tail  of  the  helicopter  in  the  positive  counterclockwise  di¬ 
rection  toward  the  nose. 

The  helicopter  is  assumed  to  be  moving  with  a  constant 
forward  velocity,  V. 

The  blade  flapping  equation  of  motion  will  be  derived  for 
a  blade  element  dm  which  is  located  at  the  distance  r  measured 
from  the  blade  hinge .  The  overall  length  of  the  blade  will  be 
refered  to  as  R. 


Flapping  Hinge  Moments 

Using  the  assumption  of  small  helicopter  roll,  pitch,  and 
yaw  rates  compared  to  the  angular  velocity  ft  of  the  rotor,  along 
with  a  negligible  offset  of  the  blade  hinge  from  the  rotor  hub 
center,  the  blade  hinge  can  be  treated  as  an  inertially  fixed 
point  in  space .  These  assumptions  allow  for  the  sum  of  the  mo¬ 
ments  about  the  hinge  to  equal  the  time  rate  of  change  of  its 
moment  of  momentum.  The  sum  of  the  moments  will  be  taken  in 


the  Z. 


direction . 


The  external  forces  on  a  blade  element  are  depicted  in 
Figure  2  and  consist  of  the  element  lift  force,  weight,  and 
centrifugal  force  generated  by  the  angular  velocity  of  the 
blade.  Aerodynamic  drag  on  the  element  is  not  included  since 
it  does  not  contribute  to  the  flapping  hinge  moment. 

Using  the  assumption  of  a  small  angle,  the  inertial 
force  about  the  flapping  hinge  moment  is; 

“Singe  =  "“■£  (1> 

which  acts  at  a  distance  r  from  the  hinge.  The  centifugal 
force  on  the  blade  element  is; 

mr  a2  (2) 

This  force  component  acts  with  a  moment  arm,  r/3  ,  on  the 
hinge.  Gravitational  force  of  magnitude,  g  dm,  acts  at  a  dis¬ 
tance  r  from  the  hinge. 

The  differential  thrust  generated  by  the  blade  element  is 
expressed  as; 

dT  =  Ldr  (3) 

where  L  represents  the  lift  generated  by  the  element.  This 
force  acts  with  a  moment  arm,  r,  on  the  blade  hinge. 

Summing  the  external  force  moments  to  the  time  rate  of 
change  of  the  flap  hinge's  moment  of  momentum  (Equation  (1)), 
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yields  the  expression; 


fR mr  £  rdr  +  /**m fl2r(r  £  )dr  + /*R  rdT  +  fR  rgdm  =  0 

Jq  */0  •'O  •'O 


(4) 


Recognizing  that  the  blade  moment  of  inertia  about  the  Zh^nge 
axis  as ; 


I  =  fR  mr2dr  (5) 

JQ 

Equation  (4)  can  be  rewritten  as; 

I  ♦  u/i){r  rdT  -  J  rgdmj-  (6) 

This  is  the  initial  form  of  the  blade  flap  angle  equation  of 
motion.  Further  refinements  are  necessary  to  make  Equation 
(6)  more  useful. 

Blade  Aerodynamics 

The  aerodynamic  force  on  the  blade  element  will  next  be 
examined  using  the  quantities  as  depicted  in  Figure  3* 

The  forward  velocity  on  the  element  is  depicted  as  U. 

This  velocity  has  a  component  in  the  direction  u^.  It 

also  has  a  radial  component,  uR,  in  the  X^a(je  direction. 

The  blade  is  at  an  ane-le  of  attack,  a  ,  to  the  velocity 
vector.  This  angle  is  the  sum  of  the  blade  pitch  angle  to  the 
horizontal  X,  .  ,  Y,  .  dane,  8  ,  and  the  velocity  vector 

angle  to  the  Xh^nge,  Yh^n^.g  plane,  It  is  mathematically; 


Assuming  small  <£,  the  following  trigonometric  relationships  can 
be  applied; 

sin  <j>  =  tan  <£  =  Up/uT  (8) 

where  uT  represents  the  tangential  velocity  to  the  blade  element 
Using  Equation  (8),  Equation  (7)  is  rewritten  as; 

a  =  e  +  Up/uT  (9) 

The  lift  equation  for  the  blade  element  is  written  as; 

L  =  ip  U2cC1  (10) 

where  c  represents  the  chord  of  the  airfoil  section  and  C1 
represents  the  lift  coefficient  for  a  given  a  .  C-^  can  be 
approximated  by; 

C1  =  a  a  ( 11 ) 

where  a  represents  the  versus  a  curve  slope.  Using  Equation 
(9)f  Equation  (11)  is  written  as; 


C1  =  a  {«  +  Up/uT  } 


(12) 


By  recognizing  that  for  small  4>  ,  U=uT  and  the  application  of 
Equation  (12),  Equation  (10)  is  transformed  into? 


L 


=  i p  u^ca  { 8 


(13) 


The  differential  thrust,  as  expressed  in  Equation  (3).  using 
Equation  (13)*  is  now  written  as; 


dT  =  i p  u^ca  +  Up/uT}  dr 


(14) 


The  terms,  uT»  up,  and  of  Equation  (14)  will  next  be  further 
refined . 

As  previously  mentioned,  uT  represents  the  velocity  in  the 
Yhinge  direction,  up  represents  the  velocity  component  in  the 
Zhinge  direc’ti°n*  and  ur  represents  the  velocity  in  the  x^^ade 
direction. 

The  term,  u,p,  is  composed  of  two  forward  velocity  components 
of  the  blade  through  a  stationary  air  mass.  These  components  of 
uT  are  the  tangent  velocity  of  the  blade  element,  r SI,  due  to 
the  angular  rotation  of  the  blade,  and  the  tangent  component  of 
the  helicopter's  forward  velocity  to  the  blade,  Vsin  \f>  .  Thus, 
uT  can  be  mathematically  expressed  as; 


u,p  =  r  £1  +  Vsintj/ 


(15) 
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The  radial  velocity  component,  uR,  along  the  Xbiad6  axis 
is  composed  of  the  normal  velocity  component  of  forward  heli¬ 
copter  motion  and  is  expressed  as; 

uR  =  Vcos^  (16) 

The  perpendicular  velocity  component  of  the  blade,  up,  is 
composed  of  the  velocity  component  in  the  g  direction  due 

to  the  flap  angle  velocity,  f3  ,  and  the  induced  velocity  due  to 
downwash,  v,  as  predicted  by  finite  airfoil  theory.  By  in- 
volking  the  assumption  of  £  being  small,  up  can  be  expressed  as; 

up  =  v  -  /3  Vcos  f  -  (17) 

Nondimensionalizing  the  velocity  components  by  the  term, 

Ril  ,  Equations  (15)-(17)  can  be  written  as; 


UT  =  x  +  ^.sin«|/ 

UR  =  fiCQSrj/ 

Up  =  X  -  ft /3  cos  \\t  -  x/3/Q. 


(18a) 

(18b) 

(18c) 


where  X  represents  the  inflow  ratio,  ft  is  defined  as  the  blade 
advance  ratio,  and  x  represents  the  nondimens ional  distance 
along  the  blade.  The  terms,  ft  and  x  are  represented  by; 


fi  =  V/(Rfl  ) 

(19) 

X 

it 

§ 

(20) 

Nondimens ionalizing  Equation  (14)  by  the  terms,  (|/»c a) 
and  R,  along  with  substituting  the  nondimens ional  velocity 
components  of  Equations  (I8a)-(l8c),  yields  the  nondimens ional 
aerodynamic  moment  for  the  blade; 

[ipcafl2R4]/u^x(Up/UT  +  8  )dx  (21) 

For  this  study,  the  local  pitch  of  the  blade  was  considered 
to  be  composed  of  four  components.  The  first  component,  Qq, 
represents  the  pitch  angle  of  the  blade  root.  The  second  com¬ 
ponent,  0-^,  represents  a  built  in  linear  twist  angle  along  the 
blade's  longitudinal  axis.  The  third  and  fourth  components,  9 

s 

and  B  c,  represent  the  cyclic  pitch  control  terms.  These  control 
terms  are  a  function  of  blade  azimuth  angle  and  represent  the 
change  in  blade  pitch  necessary  for  directional  control  of  the 
thrust  vector  across  the  rotor  plane  for  helicopter  horizontal 
and  vertical  control. 

The  pitch  angle,  6  ,  can  thus  be  represented  as; 

9  =  9  Q  +  x01  +  flgsin\|f  +  0ccos  \|f  (22) 

Substitution  Equations  (22)  and  (18c)  into  Equation  (21) 
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To  establish  the  integration  boundaries  of  Equation  (23). 
two  concepts  must  be  addressed. 

The  first  concept  is  the  blade  tip  loss  factor,  B.  Due 
to  the  limitations  of  blade  element  theory  in  predicting  the 
performance  of  a  blade  of  finite  length,  adjustments  must  be 
made  to  compensate  for  the  loss  of  lift  generating  ability  at 
the  blade  tip.  This  adjustment  is  made  by  defining  the  blade 
tip  loss  factor  as  a  limit  of  integration  along  the  nondimen- 
sional  blade  longitudinal  direction.  The  ad  hoc  selected  nu¬ 
merical  value  for  B  is  0.97,  which  gives  good  correlation  with 
experimental  data  (28j60). 

The  second  concept  which  must  be  addressed  is  the  defi¬ 
nition  of  flow  regions  across  the  rotor  plane  (28» 152,41 , 8) . 
These  regions,  as  depicted  in  Figure  4,  are  refered  to  as  region 
1  for  normal  flow,  region  2  for  mixed  flow,  and  region  3  for 
reversed  flow  (38) . 

The  first  region  is  considered  normal  flow  because  the 
air  approaches  the  blade  element  from  the  leading  edge.  This 
region  is  defined  from  where  the  azimuth  angle  equals  0  degrees 
to  where  it  equals  180  degrees.  The  third  region  is  called 
reversed  flow  because  the  air  approaches  the  blade  from  the 
trailing  edge.  The  angular  boundaries  for  this  region  are  de¬ 
fined  as; 


^  =  7r+«  to  f  =  2  ir  -  « 


(24) 


where  «  is  defined  as; 


In  the  reversed  flow  region,  the  lift  coefficient  is  assumed 
to  be  the  negative  of  the  lift  coefficient  for  the  given  angle 
of  attack  in  normal  flow  (4l). 

The  second  region  of  flow  across  the  rotor  is  called  mixed 
flow  because  it  represents  a  combination  of  normal  flow  across 
the  tip,  and  reverse  flow  close  to  the  hub.  The  angular  regions 
that  define  its  boundaries  ares 

f  =  ir  to  i r  +  t  and  if/  -  2  w  -  «  to  2  ir  (26) 

To  handle  this  condition,  the  appropriate  aerodynamic  expres¬ 
sions  must  be  evaluated  for  each  region.  This  means  that  the 
rotor  flow  model  must  be  an  integral  part  of  the  blade  flapping 
simulation  model.  There  have  been  alternative  approaches  to 
this  condition  (11)  which  gives  approximately  the  same  azimuth 
angle  boundaries  for  the  different  flow  conditions. 

Other  flow  effects  which  can  be  considered,  but  are  out 
of  the  scope  in  this  study,  are  those  due  to  transonic  effects 
(23),  wake  effects  (25,26),  and  atmosopheric  turbulence  (15,37). 

For  the  mixed  flow  region,  the  normal  velocity  approaching 
the  trailing  edge  of  the  blade  is  equal  to  the  normal  ve¬ 
locity  approaching  the  leading  edge  of  the  blade  at  the  non- 
dimensional  blade  location,  -  fi  sin  This  establishes  the 
boundary  for  normal  and  reversed  flow  on  the  blade  when  it  is 


in  the  mixed  flow  region. 


Using  the  definitions  of  tip  loss  factor  along  with  the 
three  distinct  flow  regions,  the  nondimens ional  aerodynamic 
moment  equation  has  for  its  integration  boundaries; 


Region  It  f 
Jn 


Region  Zt  J 


sin^ 


^.sm  \)r 


Region  3»  -  fB 
J0 

Using  the  integration  boundaries  along  with  combining  terms 
in  /3  and  /3  ,  Equation  (23)  can  be  expressed  as; 

{Ef(  *>  -  *k(  *>  -  (£/a)c(  *)}  (27) 


where 


K(  tff )  -  ft  J  U^xcos^dx 
C(  *)  -  /uTx2dx 


Equation  (28)  represents  the  aerodynamic  spring  charac¬ 
teristics  on  the  flap  moment.  Equation  (29)  represents  the 
aerodynamic  damper  characteristics  of  the  blade  on  the  flap 
moment.  The  term,  f(^),  represents  a  forcing  function  due  to 
aerodynamic  forces  resulting  from  applied  pitch  controls,  initial 


blade  pitch,  built  in  linear  twist,  and  inflow  effects.  Thus, 
f  (  «//  )  represents  all  the  aerodynamic  and  control  effects  built 
around  the  flapping  dynamics  of  the  blade . 


The  expression  for  K(  ^  )  and  C(  \fr  )  have  been  determined 
for  each  flow  region  (33.36,21),  and  are  presented  in  Table  I. 

Substituting  Equation  (27)  into  Equation  (6)  yields  the 
flapping  equation  of  motions 

&  +  CLZ0  =  [~ 2T~~Ca]  {lf(  *  )  “ 

-  (4  /a  )c(  *  )} 

Defining  the  blade  lock  number  as? 


-  /3K  (  *  ) 


rn 

/  rgdm 

Jc\ 


(30) 


and  introducing  a  flap  restraint  spring  with  the  spring  constant 
Kq,  which  produces  a  moment  at  the  flap  hinge  of; 


Ko* 


(32) 


to  counteract  the  aerodynamic  moment  at  the  flap  hinge,  Equation 
(30)  can  be  written  as; 


0  +  C(  *  )  $  +  |  P2  +  K(  *  )  ]  £  = 

-r2^~  +  /rgdm  (33) 
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MOMENT  FUNCTIONS 


where ; 


,  K 

P  =  1  +-~  (34) 

a  i 

Equation  (34)  represents  the  natural  frequency  of  the 
flapping  motion  when  the  helicopter  is  hovering.  Equation  (31) 
represents  the  ratio  of  the  aerodynamic  forces  to  inertial  forces. 

By  choosing  an  appropriate  hinge  spring  constant  Kq,  a 
hingeless  rotor  blade  can  also  be  examined  through  Equation  (33)- 

For  this  stability  study,  only  the  homogeneous  part  of  Equation 
(33)  was  considered.  The  only  situation  where  the  forcing 
functions  of  the  heterogeneous  part  of  Equation  (33)  would 
impact  stability,  would  be  when  the  homogeneous  part  was 
unstable  (21).  For  a  stable  system,  the  effects  of  the  forcing 
function  terms  would  be  a  different  equilibrium  pitch  and  flap 
angle  position  other  than  zero  at  some  azimuth  angle. 

Stability  Model 

To  eliminate  the  dependence  of  rotor  angular  velocity, 

SI,  on  the  solution,  a  coordinate  transformation  was  made. 

Assuming  a  constant  rotational  velocity  and  using  the  chain 

*  «• 

rule,  /9  and  $  were  transformed  to  be  a  function  of  azimuth  angle 
rather  than  time.  This  transformation  took  the  form; 


dt  d  tfr  ^ 


(35) 


and; 
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r-1. p.  p 


it  [If]  -  °2  t^r2  =  *V 


(36) 


The  new  form  of  the  flapping  equation  after  making  the  sub¬ 
stitution  of  Equations  (35)  and  (36),  is; 


0"  +  |’c(*)0'+  |  [-|  P2  +  K(  *  )]  0  =0 

Converting  Equation  (37)  into  state  space  representation 
yields  the  expression; 


1  [? p2  +  K<*>] 


(37) 


V 

/. 

1 - 

— — • 

oa.. 

(38) 


Equation  (38)  was  used  as  the  simulation  model  for  the 
blade  flapping  motion.  The  rate  of  change  of  flap  angle,  0  , 
was  initially  set  at  zero,  with  the  flap  angle  perturbed  from 
the  equilibrium  position  at  the  zero  azimuth  degree  location. 

Had  the  matrix  coefficients  of  Equation  (38)  been  constant 
values,  standard  classical  and  modern  control  procedures,  such 
as  root  locus  and  Bode  analyses,  could  have  been  used  for  sta¬ 
bility  determination  and  control  design.  However,  with  the 
periodic  coefficients,  these  techniques  could  not  be  used  and 
a  different  approach  had  to  be  taken.  This  approach  centers 
around  Floquet  theory  which  will  be  addressed  in  the  next 
section . 


25 


III.  Floquet  Theory 


From  Section  II,  it  was  demonstrated  that  the  helicopter 
blade  flapping  equation  of  motion  represented  a  nonlinear, 
time-periodic  differential  equation.  It  was  further  demon¬ 
strated  that,  through  the  use  of  appropriate  assumptions  and 
the  physical  characteristics  of  the  blade,  these  equations  can 
be  transformed  into  linear,  time-periodic  differential  equations 
such  as  Equation  (38).  For  low  speed  forward  flight  where  the 
advance  ratio  is  relatively  low,  it  can  be  assumed  that  the 
coefficients  of  the  flapping  equation  of  motion  are  constant 
values.  However,  for  the  faster  forward  velocity  case,  the 
periodic  nature  of  the  differential  equation  returns. 

The  most  straightforward  method  of  dealing  with  the 
stability  of  this  type  of  periodic  system  is  through  the  appli¬ 
cation  of  Floquet  theory  (9.21,36).  The  major  utility  of 
Floquet  theory  is  to  determine  the  stability  of  the  dynamic 
system.  For  the  helicopter  blade  application,  the  theory  has 
been  used  only  for  this  purpose  and  not  as  a  design  synthesis 
tool  for  feedback  control. 

The  purpose  of  this  section  is  to  lay  out  the  general 
foundations  of  Floquet  theory.  The  consequences  of  the  theory 
will  be  addressed  along  with  some  insights  on  how  the  theory 
can  be  applied  for  feedback  control  synthesis. 

Foundations  of  Floquet  Theory 

To  begin  the  development  of  the  basic  ideas  of  Floquet 
theory,  a  periodic  system  will  be  defined  as; 


where  the  matrix  A(t  )  is  periodic  over  some  time  interval 
T.  The  solution  to  this  differential  equation  is; 

X(t)  =  $(t,0)  X(0)  (40) 

where  3>(t,0)  represents  the  state  transition  matrix.  The 
solution  to  the  state  transition  matrix  is  obtained  by  solving 
the  differential  equation; 

ft  4>(t,0)  =  Aft)  *(t,0)  (41) 

which  has  the  identity  matrix  for  its  initial  condition. 

Floquet  theory  states  (4,9)  that  the  state  transition 
matrix  can  be  written  as; 

<J>(t.O)  =  F(t)  exp{[j]  t}F_1(0)  (42) 

where  t  is  any  time  greater  than  zero.  In  Equation  (42), 
the  matrix  [j]  represents  a  constant  matrix  which,  in  most 
cases,  is  constructed  in  Jordan  normal  form.  The  diagonal 
entries  of  j^jj  are  refered  to  as  the  Poincare'  exponents, 
which  are  similar  to  the  eigenvalues  of  a  constant  coefficient 
differential  equation.  The  matrix  F(t)  is  periodic  with  the 
same  period  of  the  dynamic  system. 


It  should  be  noted  here  that  Floquet  theory  does  not 
yield  the  solution  to  the  problem,  but  gives  information  about 
the  form  and  properties  of  the  solution  (9).  The  knowledge 
over  one  period  of  time  determines  the  solution  to  the  homo¬ 
geneous  system  of  Equation  (39)  through  Equation  (42).  Thus, 
solving  the  Floquet  problem  for  all  time  requires  finding  the 
constant  matrix  [jj  and  the  periodic  matrix  F(t)  over  one 
period . 

Since  F(t)  is  periodic,  the  elements  of  F(t)  at  the  be¬ 
ginning  of  the  period  equals  the  elements  at  the  termination 
of  the  period.  Using  this  information  to  set  the  initial 
conditions,  solving  Equation  (42)  yields  the  solution; 


<X>  (T,  0)  =  F( 0)  exp{[j]  t}  F“1(0) 


(43) 


In  this  solution,  <X>(T,0)  is  refered  to  as  the  monodromy  matrix. 
The  eigenvectors  of  this  monodromy  matrix  represents  the  column 
vectors  of  the  periodic  matrix  F(0)  or  F(T).  Since  4(T,0) 
and  expjjjj  ^  are  similar  matrices,  the  eigenvalues  of  the 
monodromy  matrix  4>(T,0)  are  related  to  the  Poincare'  exponents 
through  the  relationship; 


Xi  =  exp(  w.T) 


As  previously  mentioned,  Floquet  theory  just  gives  infor¬ 
mation  about  the  form  and  properties  of  the  solution,  and  not 
the  solution  itself.  In  order  to  obtain  the  solution,  the 


differential  equation  for  the  state  transition  matrix,  Equation 
(4l),  must  be  directly  integrated  using  numerical  methods. 

This  integration  involves  simulating  the  state  transition 
matrix  for  one  period  starting  with  the  appropriate  initial 
conditions  of  the  identity  matrix.  The  topic  of  numerical 
methods  used  for  this  study  is  presented  in  Appendix  A. 

From  a  historical  point  of  view,  it  is  interesting  to 
mention  the  development  of  helicopter  blade  stability  analyses. 
Prior  to  the  introduction  of  the  digital  computer,  the  heli¬ 
copter  blade  flapping  motion  was  modelled  using  Hill's  equation 
(21).  The  solutions  to  this  class  of  differential  equations 
were  precomputed  and  tabulated  assuming  a  standard  form.  The 
blade  flapping  motion  equations  were  converted  to  fit  this 
general  form,  and  stability  was  determined  through  the  use  of 
Floquet  theory.  This  approach  generally  did  not  allow  for 
many  parametric  studies  to  be  conducted  and  in  general,  the 
models  themselves  had  to  be  less  complicated  allowing  for 
mathematical  simplicity.  When  the  general  equations  including 
the  aeroelastic  effects  on  the  blade  were  developed  (22),  the 
digital  computer  was  used  to  simulate  the  motion  of  the  blade 
in  forward  motion  (23).  Stability  was  determined  by  simulating 
the  motion  of  the  blade  for  several  revolutions  until  the 
growth  of  the  oscillations  of  the  blade  could  be  determined. 
Later,  this  technique  was  replaced  by  Floquet  theory. 

The  eigenvectors  of  the  raonodromy  matrix  are  calculated 
as  follows.  The  solution  is  obtained  by  substituting  Equation 
(43)  into  Equation  (4l),  which  yields? 


55  [?(t)eJtF"1(0)}-  =  A(t)F(t)eJtF-:l(0) 


(45) 


Carrying  through  the  matrix  differentiation  of  Equation  (45) 
results  in  the  differential  equation; 

55  F(t)  =  A(t)F(t)  -  F(t)[  j]  (46) 

The  initial  conditions  used  for  the  integration  of  Equation 
(46)  are  the  eigenvectors  of  the  monodromy  matrix.  After  the 
integration  is  completed  for  one  period,  the  individual  values 
of  the  elements  of  matrix  F(t)  as  a  function  of  time  are  con¬ 
verted  into  an  appropriate  Fourier  series . 

The  stability  of  the  dynamical  system  can  be  determined 
through  Floquet  theory  by  either  looking  at  the  magnitude  of 
the  characteristic  multipliers  (eigenvalues  of  the  monodromy 
matrix),  or  the  algebraic  sign  of  the  associated  Poincare' 
exponents.  The  conditions  for  dynamic  stability  of  the  system 
are  thus ; 

Vi  Xi(real))2  +  (  X  i(  imaginary)  )2  <1  (47) 


or; 


<ui  <  0 


(48) 
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For  a  feedback  controller,  the  system  model  is  described 
as; 

jj.  X(t)  =  A(t)X(t)  +  B(t)u(t)  (49) 

Making  the  transformation  into  Floquet  modal  form  (44); 

X(t)  =  F(t)  rj  (50) 

ft  x(t)  ■  {ft  F<t)}  7  *  F(t)  {4^,}  (51) 

and  substituting  Equation  (46)  into  Equation  (51)  yields; 

X(t)  =  {A(t)F(t)  -  F(t)j}  v  +  F(t)^r,  (52) 

Substituting  Equations  (50)  and  (52)  into  Equation  (49),  the 
system  model  in  modal  form  becomes; 

-  JV  +  F"1(t)B(t)u(t)  (53) 

For  a  two  dimensional  case  in  phase  variable  state  space 
form,  as  in  Equation  (38),  the  control  matrix,  B(t),  is  repre 
sented  by; 


Had  the  system  model  been  constructed  in  control  canonical 
form,  the  control  matrix  B(t)  could  have  been  defined  as; 

(1  0)T  or  (1  1)T 

However,  in  phase  variable  form,  the  only  definition  allowed 
was  that  of  Equation  (5*0,  which  represents  a  torque. 

The  state  feedback  term,  u(t),  is  represented  by; 

u(t)  =  (k^t)  .  kg( t)  (55 a) 

=  K(t),  (55b) 

and; 

G(t)  =  F”1 ( t )B  (56a) 

=  (gx(t)  g2(t))T  (56b) 

The  matrix  G(t)  of  equation  (56a)  represents  the  internally 
generated  mode  controllability  matrix  of  the  system  model  as 
a  result  of  the  product  of  the  inverse  eigenvector  matrix, 
F~1(t),  and  the  control  matrix  B. 

In  Equation  (556),  K(t)  represents  the  gain  variables 
which  can  be  selected  by  the  control  designer.  Thus,  through 
K(t),  the  control  designer  has  the  ability  of  single  or  double 


modal  control. 


Single  Mode .  Assuming  that  the  second  modal  variable  is 
of  interest  to  control,  K(t)  takes  on  the  form; 

[o  k2(t)]  (57) 


Substituting  this  matrix  into  Equation  (53)*  along  with  the 
G(t)  matrix  as  described  in  Equation  (56b),  yields  the  system 
model  in  modal  form  as; 


d_ 

dt 


V 


k2(t)g1( t) 


«2  +  k2(t)g2(t) 


(58) 


In  the  two  dimensional  case  for  scalar  control,  the 
selection  of  the  K(t)  matrix  will  either  make  the  matrix  of 
Equation  (58)  in  either  upper  triangular  form,  such  as  it 
currently  is,  or  in  lower  triangular  form.  Approaching  Equation 
(58)  as  a  Floquet  type  problem,  the  new  Poincare'  exponents 
of  the  system  will  be; 


ul'  =  wl  (59) 

«2’  =  «2  +  (jr) f  k2(t)g2(t)dt  (60) 

Thus,  for  this  particular  case,  the  first  Poincare'  ex¬ 


ponent  remains  unchanged  while  the  second  Poincare  ’  exponent 
is  shifted  by  the  integral; 


^  rTk2(t)g2(t)dt 


(61) 


The  upper  or  lower  triangular  form  of  Equation  (58)  will 
not  prevail  for  higher  dimensional  systems?  however,  the 
sparcity  property  of  the  resulting  matrix  of  these  systems  for 
single  mode  control  will  make  the  properties  as  described  by 
Equations  (59)  and  (60)  applicable.  An  example  of  this  cam  be 
seen  in  a  third  order  system  where  the  K(t)  matrix  is  selected 
as  j 


K(t)  =  [o  k2(t)  o] 


The  modal  system  differential  equation  would  be? 


(62) 


coj  g1(t)k2(t) 

0  "2  +  g2^k2^ 

0  g^(t)k2(t) 


The  new  Poincare'  exponent  «2'  is  given  by; 


iiir 


w2  +  (T}/o  *2(t)k2(t)dt 


(63) 


(64) 


and  the  remaining  exponents  are  unchanged. 

Again,  there  is  only  a  change  in  one  Poincare’  exponent 
while  the  other  exponents  remain  unchanged.  However,  due  to 
the  cross  coupling  terms,  there  would  be  a  change  in  the  system 
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e igenvectors . 


The  solution  of  Equation  (58)  is  expressed  as; 


•7  ( t )  =  7(0) exp  J 


<u1  k2(t)g1(t) 


0  «2  +  k2(t)g2(t) 


Single  Scalar  Modal  Control .  The  first  case  to  be  examined 
is  where  the  fourier  series  transf ormation  of  the  column  matrix 
components  resulting  from  Equation  (58a)  can  be  separated  into 
a  constant  D.C.  component  and  a  periodic  A.C.  component.  For 
the  present  example,  this  results  into  the  condition; 


gl(t)  =  «1  D.C.  +  gl(t)A.C. 
g2(t)  =  g2  D.C.  +  g2(t)A.C. 


(66) 

(67) 


Using  this  expression  for  g^(t)  and  g2(t),  along  with  a  constant 
gain  of  k2  in  Equation  (57).  yields  the  expression; 


dt 


(68) 
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The  stability  of  17  (t)  is  determined  by  the  eigenvalues 
of  the  matrix? 


W1  k2gl  D.C. 

0  "2  +  k2g2  D.C. 


(69) 


These  by  inspection  are; 

W1  “2  +  k2g2  D.C. 

Thus,  for  scalar  control,  the  new  Poincare'  exponent  lo¬ 
cation  can  be  achieved  through  the  expression; 


w , '  =  w  o  +  k-g 


2&2  D.C. 


(70) 


For  the  case  where  g2(t)  does  not  have  a  constant  part  in 
its  Fourier  series,  this  method  is  still  applicable.  For  this 
case,  k2(t)  is  not  considered  a  constant  value,  but  rather  can 
take  on  the  form,  such  as  (*0; 


k2(t)  =  k  sin  2 -3  - 


(71 ) 


Through  the  multiplication  of  k2(t)g2(t),  the  product  will 
generate  a  constant  term  through  standard  trigonometric  i- 
dentities .  This  term  would  be  kc  ,  where  c  represents  the 

rig 
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un 


n  sine  term  of  the  Fourier  series  expansion  of  g2(t) 
Thus,  the  new  Poincare'  exponent  would  be  given  as; 


kc 


n. 


w2’  =  "2  +  T 


(72) 


This  procedure  would  also  work  if  a  cosine  function  was 

used  for  Equation  (72),  where,  in  this  situation,  c  would 

n2 

represent  the  n^h  cosine  term  in  the  Fourier  series  expansion 
of  g2( t) . 

Scalar  Control  of  Two  or  More  Modes .  This  case  examines 
where  it  is  desired  to  control  more  than  one  mode.  The  approach 
is  applicable  for  both  real  or  complex  Poincare'  exponents. 

The  gain  matrix  K(t)  of  Equation  (55a)  is  represented  by 
a  row  vector  composed  of  individually  selected  gains  per  mode; 

K(t)  =  [k1(t),  k2(t) . kn(t)]  (73) 

Substituting  K(t)  along  with  Equation  (5*0  and  (56b),  into 
Equation  (53)  yields  the  result; 


d 

dt  ^ 


»1  +  g1(t)k1(t) 
g2(t)k1(t) 


g1(t)k2(t) 


(02  +  g2(  t)k2(t) 


(74) 


This  equation  is  applicable  if  k1(t)  and  k2(t)  represent 
constant  gains,  or  gains  which  can  be  expressed  as  a  sine  or 


r 


0 


r 


l: 


cosine  function  of  the  mode  and  the  period.  With  the  gains 
g1(t)  and  g2(t)  represented  as  a  Fourier  series,  Equation  (74) 
can  be  written  as? 


v  =  A.1  17  +  A2(t)i7  (75) 

where  A1  represents  a  constant  matrix  and  A2(t)  represents  a 
matrix  with  periodic  terms . 

The  shortcoming  of  this  approach  is  that,  by  choosing  A1 
for  negative  eigenvalues,  this  selection  does  not  guarantee  a 
stable  system  (4,42).  There  is  no  general  way  of  determining 
gains  such  that  Equation  (75)  has  predictable  Poincare’  ex¬ 
ponents.  With  this  limitation,  the  only  alternative  to  control 
design  is  to  go  through  the  previously  described  design 
synthesis  of  numerically  solving  the  Floquet  system  to  determine 
the  proper  gains.  However,  there  is  some  insights  for  choice 
of  gains  to  be  offered  in  this  approach. 

Using  the  expression  of  a  matrix  derivative  from  linear 
algebra  for  a  first  order  ordinary  differential  equation; 

—  D(t)  =  trace  |A(t) j  D(t)  (76) 

where  D(t)  represents  the  determinant  of  the  monodromy  mccrix, 
and  A(t)  represents  the  matrix  of  Equation  (74).  By  evalu¬ 
ating  the  expression  for  time  equalling  one  period  and  using 
Poincare'  exponents,  the  solution  to  Equation  (76)  can  be 


r: 


»» 


••  ,*  *  .*  *  »  *  _»*•' 
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expressed  asj 


z 


a> 


i 


■jr  f  trace  A(t)  dt 
1  J0 


(77) 


For  the  two  dimensional  modal  case,  this  expression  becomes! 


wl’  +  “2*  "  W1  +  w2 

+  "T  /0  {kl(t)gl(t>  +  k2(t)g2(t)}  dt  (?8) 

Assuming  that  g^(t)  and  g2(t)  have  constant  terms  in  their 
Fourier  series,  constant  gains  can  be  used  and  Equation  (78) 
can  be  written  as; 

wl‘  +  w2’  =  W1  +  w2  +  klgl  +  k2g2  (79) 

This  expression  gives  the  sum  of  the  new  Poincare '  ex¬ 
ponents  but  it  cannot  predict  what  the  individual  values  of 
the  exponents  would  be.  A  systematic  design  approach  thus 
is  required  to  determine  the  individual  gains  k^  and  k 2  to 
give  the  desired  system  performance.  It  should  be  mentioned 
that  this  approach  can  be  used  for  the  case  of  complex  conju¬ 
gate  exponents  to  predict  the  real  part  of  the  conjugates. 
Recapitulation 

For  the  solution  to  the  control  of  a  time  periodic  dynamic 


system,  the  following  steps  are  taken. 


First,  the  original  Floquet  system  must  be  solved.  The 
first  step  to  the  solution  is  to  solve  the  monodromy  matrix 
by  integrating  the  dynamic  system  of  equations  over  one  period. 
From  the  monodromy  matrix,  the  Poincare’  exponents  and  the 
stability  of  the  system  can  be  determined.  The  solution  of 
the  monodromy  matrix  also  yields  the  initial  values  for  the 
time  dependent  eigenvectors  of  the  periodic  system. 

Next,  the  time  dependent  eigenvectors  are  determined  from 
direct  integration  of  Equation  (46).  This  step  begins  with 
using  the  eigenvectors  of  the  monodromy  matrix  as  the  initial 
condition  to  Equation  (46).  Equation  (46)  is  then  simulated 
over  one  period  with  the  eigenvectors  evaluated  over  equally 
spaced  time  intervals.  Finally,  these  values  are  converted  into 
a  Fourier  series  representation  of  the  eigenvectors  as  a  function 
of  time,  or  as  in  this  study's  case,  a  function  of  azimuth 
angle . 

When  going  through  a  control  synthesis  of  the  time 
periodic  dynamic  system,  a  control  philosophy  must  first  be 
selected.  The  choices  of  control  as  presented  in  this  study 
are  scalar  control. 

The  approach  to  the  control  design  process  is  to  first 
form  the  appropriate  gain  expressions.  This  is  done  by  forming 
the  mode  controllability  matrix  through  Equation  (56a),  using 
the  information  from  the  eigenvectors  and  the  desired  com¬ 
position  of  the  B  matrix.  Next,  the  proper  k  gain  values  are 
selected  using  the  information  from  the  monodromy  matrix  and 
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the  desired  pole  locations.  Finally,  the  resulting  feedback 
control  law  is  placed  into  the  simulation  model  to  verify  the 
achievement  of  the  desired  pole  placement. 

Appendix  A  addresses  the  variety  of  numerical  methods 
and  computer  software  used  for  the  individual  steps  in  this 
design  synthesis. 


IV. 


ults 

This  section  will  present  the  study  results  of  the 
helicopter  blade  flapping  problem  using  as  the  mathematical 
model,  Equation  (38). 

The  case  for  no  feedback  control  will  first  be  explored. 
Next,  a  feedback  controller  for  state  feedback  will  be  ex¬ 
amined.  In  this  case,  Y =8  which  is  the  lock  number  for  a 

p 

full  size  blade  (18).  P  was  chosen  as  1.0  which  represents 
an  unconstrained  blade.  This  turns  out  to  be  the  worst  case 
condition  for  blade  stability.  Other  values  used  for  this 
analysis  were  /*  =  2.4,  and  B=0.97* 

No  Feedback  Case 

The  models  were  initially  run  without  feedback  for  com¬ 
parisons  to  other  results.  For  these  cases,  Y  was  varied  from 

2 

2  to  16  by  increments  of  2.  P  was  varied  from  1.0  to  1.8  by 
step  sizes  of  0.2  for  each  Y  .  In  all,  this  represented  40 
cases.  The  Poincare'  exponents  were  derived  for  each  case 
using  the  computer  models  as  described  in  Appendix  A.  The 
results  of  these  runs  are  presented  in  Table  II.  The  shaded 
region  of  Table  IT  represents  the  cases  where  the  system  is 
unstable.  Examples  of  a  stable  and  unstable  response  to  a  step 
input  of  one  degree  is  presented  in  Figure  5- 

A  comparison  of  the  stability  boundary  predicted  by  these 
cases  and  those  of  previous  studies  (41,36),  can  be  found  in 
Figure  6.  As  can  be  seen  in  this  figure,  there  is  close 
comparisons  between  the  results  of  this  study  to  other  studies. 
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TABLE  II 

40  CASES 


Figure  5.  Stable  and  Unstable  Response 


Figure  6.  STABILITY  BOUNDARY 


Feedback  Case 


As  previously  mentioned,  the  condition  of  Y =  8  oper¬ 
ating  with  no  blade  flapping  restraint  was  considered  for  the 
feedback  case. 

For  this  condition,  a  root  locus  plot  of  the  Poincare’ 
exponents  with  increasing  advance  ratio  is  presented  in  Figure 
7.  From  this  plot,  it  is  easily  seen  that  for  low  advance 
ratios  which  correspond  to  low  forward  speeds,  the  roots  are 
complex  conjugate  pairs .  This  compares  well  with  the  con¬ 
ventional  control  design  approaches  where  the  flapping  motion 
dynamics  are  modelled  as  second  order  systems  with  complex 
roots.  As  the  advance  ratio  is  increased,  the  roots  migrate 
to  the  real  axis,  then  split  with  one  going  to  the  left  (more 
stable)  and  the  other  travelling  to  the  right  (less  stable). 

This  is  why,  for  the  conventional  approach  to  rotor  control, 
gain  scheduling  is  done  for  the  low  speed  case,  then  readjusted 
for  higher  speeds . 

The  time  dependent  eigenvectors  for  this  case  at  increments 
of  ten  degrees  are  presented  in  Appendix  B. 

The  Fourier  series  coefficients  of  the  eigenvector  for 
this  case  and  cases  for  Y  =  2,4,6,8,10,12,14,  and  1 6 ,  along 
with  their  Poincare'  exponents  and  the  eigenvalues  of  the 
monodromy  matrix  can  be  found  in  Appendix  C. 

The  Fourier  coefficients  for  the  elements  of  the  mode 
controllability  matrix  for  these  cases  derived  by  the  computer 
software  of  Appendix  I  can  be  found  in  Appendix  D. 
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For  7=8,  the  second  Poincare*  exponent  was  changed  from 
its  present  unstable  positive  value  to  zero.  Using  the  constant 
Aq  terms  from  Appendix  D  and  Equation  (70),  the  required  gain 
to  accomplish  this  task  was  calculated  as; 

BT  =  (  0  1  ),  k=  0.109737 

Placing  these  values  into  a  feedback  control  simulation  (Appendix 
J),  the  resulting  Poincare'  exponents  were; 

U  =  -2.050636 
w2  =  0.000239 

As  can  be  seen  by  these  results,  desired  pole  placement  of 
w2  was  obtained  with  minor  change  to  «  1 .  The  reason  for  the 
slight  change  is  due  to  the  Runga-Kutta  integration  routine 
with  associated  roundoff  error  and  step  size. 

a  comparison  of  the  response  of  this  system  to  the  case 
of  no  feedback  can  be  found  in  Figure  8. 

By  applying  this  gain  to  the  other  cases  with  different 
lock  numbers,  the  system  still  remains  stable.  This  was  done 
by  using  the  Aq  terms  from  Appendix  D  and  the  Poincare'  ex¬ 
ponents  from  Table  II.  The  comparison  of  the  old  Poincare' 
exponents  to  the  achieved  exponents  is  found  in  Table  III. 


V .  Conclusion 


From  this  study,  it  was  shown  how  the  principles  of 
Floquet  theory  could  be  applied  to  the  helicopter  blade 
flapping  control  design  synthesis. 

By  first  deriving  the  appropriate  blade  flapping  equations 
of  motion  for  a  rigid  blade  with  no  twist,  a  series  of  com¬ 
puter  simulations  were  assembled.  Each  of  these  simulations 
addressed  a  different  aspect  in  solving  the  Floquet  problem. 

From  the  analysis  of  a  blade  with  Lock  number  of  8,  it 
was  shown  how  conveniently  a  feedback  controller  could  be 
designed  for  the  case  of  state  feedback.  The  conventional 
approach  to  the  same  problem  would  begin  the  design  process 
by  first  considering  low  forward  flight,  then  define  stability 
regions  for  different  forward  velocities.  The  gains  would 
then  have  to  be  readjusted  until  the  system  was  in  a  stable 
region.  This  step  has  been  eliminated  with  Floquet  theory 
application. 

This  has  been  the  first  step  in  the  process  of  realizing 
the  full  control  scheme  for  the  helicopter  rotor.  Further 
topics  which  need  to  be  addressed  are  the  coupled  blade 
flapping  with  lead-lag  degrees  of  freedom.  This  is  another  well 
suited  application  of  Floquet  theory  control  design  synthesis. 

Another  case  which  renders  itself  to  this  type  of  analysis 
is  the  case  where  the  blade  is  travelling  through  atmosopheric 
turbulence.  In  this  condition,  the  flow  model  as  described 
in  Section  II  would  be  replaced  with  a  stochastic  model. 


As  presented  in  Section  II  and  Section  III,  numerical 
methods  are  required  to  derive  the  solution  of  the  time 
periodic  blade  flapping  equation  of  motion.  These  methods 
are  required  due  to  the  fact  that  Floquet  theory  can  only 
give  information  about  the  form  and  properties  of  the  solution, 
but  not  the  solution  itself.  To  obtain  the  solution,  the 
dynamic  system  must  be  simulated  over  one  period  to  obtain 
the  proper  data  to  derive  the  monodromy  matrix,  the  periodic 
eigenvectors  of  the  system  as  a  Fourier  series  representation 
dependent  on  azimuth  angle,  and  the  periodic  gain  properties 
as  a  Fourier  series  representation  dependent  on  azimuth  angle. 

The  required  algorithms  to  do  this  analysis,  as  high¬ 
lighted  in  Section  III,  consist  of  an  integration  routine, 
eigenproblem  routine,  and  a  Fourier  series  coefficient  deter¬ 
mination  routine. 

The  purpose  of  this  section  is  to  present  the  software 
developed  to  handle  this  Floquet  type  problem.  The  software 
developed  can  be  found  in  Appendix  E,F,G,H, J.and  J,  and  will 
be  described  by  highlighting  the  specific  function  of  each 
program.  All  the  programs  have  been  written  in  Texas 
Instruments  Extended  BASIC,  which  is  a  dialect  of  the  BASIC 
programming  language,  and  is  closely  related  to  Microsoft 
BASIC-80. 

The  computer  used  to  develop  and  conduct  the  research 
of  this  study  was  the  Texas  Instrument  99-^-A  microcomputer. 


This  computer  uses  the  TMS9900  16-bit  microprocessor  which 
utilizes  a  16-bit  data  bus.  The  microprocessor  as  configured 
can  give  a  decimal  precision  of  13  to  14  digits  depending 
on  the  numerical  value,  and  an  exponential  decimal  range 
from  -128  to  +126.  These  capabilities  make  this  computer 
suitable  to  handle  the  numerical  accuracies  as  demanded  in 
the  eigenproblem  algorithm.  Other  required  computer  hardware 
to  conduct  this  research  consisted  of  48K-bytes  of  random 
access  memory,  and  a  disk  drive  for  program  storage  and  data 
manipulation. 

Integration  Routines 

Two  integration  routines  were  required  to  derive  the 
solution  to  Equation  (41)  and  Equation  (46).  Both  routines 
use  the  blade  flapping  motion  model  as  described  in  Equation 
(38). 

The  integration  technique  employed  in  both  routines  was 
a  fourth-order  Runga-Kutta  routine  (37*194, 40:182) .  The 
Runga-Kutta  routine  was  selected  because  of  its  numerical 
accuracy  and  self-starting  capability.  Other  higher  order 
Runga-Kutta  techniques  could  have  been  chosen;  however,  a 
tradeoff  was  made  between  accuracy  and  computer  running  time. 
Both  routines  integrated  with  respect  to  azimuth  angle  as 
dictated  by  Equation  (38),  with  the  step  size  of  integration 
adjusted  to  give  the  best  accuracy.  The  adjustment  of  inte¬ 


gration  step  size  will  be  presented  later  in  this  section. 

The  first  integration  routine  (Appendix  E),  was  used  to 


determine  the  monodromy  matrix  described  by  Equation  (41). 

The  function  of  the  program  and  its  subroutines  follow. 

The  executive  portion  of  the  program  first  initializes 
all  variables  used.  Afterwards,  it  calls  for  the  parameters 
of  tip  loss  factor,  advance  ratio,  Lock  number,  and  the  natural 
frequency  of  the  helicopter  blade  in  hover.  From  there,  it 
begins  the  fourth-order  Runga-Kutta  routine  to  simulate  the 
flapping  motion  of  the  blade  for  one  revolution.  Afterwards, 
the  numerical  values  of  the  elements  of  the  monodromy  matrix 
are  printed  out. 

Subroutine  FCNN  was  the  routine  which  represented  Equation 
(41).  Major  calls  from  this  subroutine  are  to  subroutine 
MATRIXA,  which  represents  the  helicopter  blade  model  of  Equation 
(38),  and  subroutine  MATRIXMULT,  which  is  a  standard  matrix 
multiplication  routine.  Subroutine  FCNN  was  called  from  the 
imbedded  Runga-Kutta  routine  of  the  executive  program. 

Subroutine  MATRIXA  contains  the  helicopter  blade  model. 
Imbedded  in  this  subroutine  where  the  equations  for  the  periodic 
coefficients  C(  t  )  and  K( ^  )  (lines  2022-2 035).  along  with 
the  determination  of  the  flow  region  (lines  2010-2016).  The 
actual  matrix  evaluation  of  Equation  (38)  was  conducted  in 
lines  2052-2055-  The  subroutine  returns  updated  values  to 
subroutine  FCNN. 

Subroutine  MATRIXMULT  was  a  standard  matrix  multiplication 
routine  (31:41).  The  subroutine  was  called  from  subroutine 
FCNN,  and  returns  its  values  to  the  same  subroutine. 

The  next  integration  routine  used  in  this  analysis  can 
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be  found  in  Appendix  G.  The  purpose  of  this  routine  was  to 
solve  Equation  (46).  This  program  also  uses  an  imbedded 
fourth-order  Runga-Kutta  routine . 

The  executive  portion  of  this  routine  was  structured 
similarily  to  that  of  the  last  program.  The  program  first 
initializes  all  variables.  The  inputs  used  in  the  program 
were  the  same  as  the  last  program  with  the  addition  of  the 
Poincare '  exponents  and  the  eigenvectors  of  the  monodromy 
matrix.  The  output  of  this  routine  were  the  elements  of  the 
eigenvector  matrix  for  every  ten  degrees  of  azimuth  angle. 

This  output  is  later  used  in  the  Fourier  series  coefficient 
analyses . 

Subroutine  FCNN  represents  the  right  hand  side  of  Equation 
(46).  This  subroutine  is  called  per  each  step  of  the  Runga-Kutta 
from  the  executive  program.  The  main  calls  from  this  sub¬ 
routine  were  to  subroutine  MATRIXA,  which  was  previously 
described,  subroutine  NIATRIXMULT  which  was  also  described,  and 
subroutine  MATRIXADD,  which  is  a  standard  matrix  addition  routine. 

Since  both  subroutines  MATRIXA  and  MATRIXMULT  have  been 
described,  the  details  of  these  subroutines  will  not  be 
presented  here. 

Subroutine  MATRIXADD  is  a  standard  matrix  addition  routine 
(31  ’•  40).  No  further  details  on  this  subroutine  will  be  pre¬ 
sented  . 

Subroutine  INPUTT  is  the  routine  which  handles  the  matrix 
input  as  required  by  the  executive  program. 


! £ 


X 


The  two  Fourier  series  analysis  routines  used  m  this 
research  can  be  found  in  Appendix  H  and  I.  Both  routines  use 
the  eigenvector  matrix  data  for  every  ten  degrees  of  azimuth 
angle  as  computed  by  the  previously  described  program.  The 
output  of  both  programs  are  the  first  19  Fourier  coefficients 
for  each  of  their  respected  cases. 

The  definition  of  the  Fourier  series  used  in  both  programs 
are  mathematically  described  in  the  following  equation  (30:479, 
19«3**.  40:148); 


[v=°s  ♦  vln  H3*] 


The  Fourier  series  coefficients  are  described  by  the  following 
equations; 


ak  =  t/0T  f(t> 


cos  3*,^  dt  k-0,1,2,3 . 18 


\  •  f I0  f(t)sln 


2  TT  tk 


dt  k=l ,  2 ,3 . 18 


As  can  be  seen  in  Equations  (81)  and  (82),  an  integration 
routine  is  also  required  in  the  software.  Because  of  its 
accuracy,  a  Simpson  integration  routine  (19:25)  was  imbedded 
in  both  programs.  It  is  important  to  note  that  the  discrete 
form  of  the  Fourier  series  analysis  was  not  used  due  to  its 
problem  of  distributing  the  Fourier  coefficients  along  all 
the  harmonics  of  the  Fourier  series.  The  continuous  model 


does  not  have  this  problem  which  now,  allows  for  the  elimi¬ 
nation  of  those  harmonics  with  lower  order  coefficients.  Both 
methods,  however,  will  give  the  same  value  of  the  constant 
component  of  the  Fourier  series  which  is  of  prime  interest  in 
the  Floquet  analysis. 

The  purpose  of  the  first  Fourier  analysis  program  (Appendix 
H)  is  to  derive  the  Fourier  series  representation  of  the 
eigenvector  matrix  for  the  homogeneous  Floquet  problem  as 
described  by  Equation  (42).  The  program  consist  of  two  parts, 
the  executive  routine  which  handles  the  data  input  and  output 
functions  along  with  the  final  stages  of  the  Fourier  analysis 
routine,  and  subroutine  FOURIER  which  handles  the  summation 
and  integration  of  the  data. 

As  previously  mentioned,  the  executive  routine  handles 
both  data  manipulation  tasks  and  completes  the  Fourier  series 
analysis.  Lines  570  and  580  computes  the  nineteen  a  and  b 
Fourier  coefficients  by  multiplying  the  summation  terms  from 
subroutine  FOURIER  by  the  Simpson  integration  constant? 

A  » 

3 

and  the  Fourier  series  constant;. 

_2_ 

T 

The  increment  of  azimuth  angle  used  in  this  analysis  was  ten 
degrees  with  the  period  of  the  system  set  at  3^0  degrees. 


(83) 


(84) 
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Subroutine  FOURIER  performs  the  integration  of  the  data 
as  dictated  by  Equations  (8l)  and  (82).  Lines  960  and  980 
are  the  locations  in  the  routine  where  this  function  takes 
place.  These  two  lines  are  also  the  location  where  the  Simpson 
integration  takes  place  through  the  multiplication  of  each 
Fourier  product  by  the  term  CONSTANT.  The  term  CONSTANT 
represents  the  coefficients  1  for  the  first  and  last  data 
entry,  4  for  odd  numbered  data  entries,  and  2  for  even  number 
data  entries.  The  determination  of  these  terms  is  conducted 
in  lines  912-916. 

The  purpose  of  the  second  Fourier  series  analysis  routine 
(Appendix  I),  is  to  evaluate  the  mode  controllability  matrix 
as  mathematically  described  by  Equation  (56a).  Although 
the  program  is  written  with  three  options  for  the  B  matrix, 
the  B  matrix  used  in  this  analysis  was  5 

BT  =  (  0  1  )  (85) 

The  structure  of  this  program  is  very  similar  to  that 
of  the  previously  described  program  with  the  addition  of  a 
matrix  multiplication  routine  and  a  matrix  inversion  routine. 
The  input  to  this  program  is  also  identical  to  the  last 
program.  The  Fourier  analysis  routines  are  identical. 

The  difference  in  the  executive  portion  of  the  program 
is  the  multiplication  of  the  inverted  eigenvector  matrix  with 
the  B  matrix.  The  product  of  this  multiplication  is  then  fed 
into  the  Fourier  analysis  routine.  This  is  accomplished  in 


the  program  at  line  locations  442-449  and  492-499.  The  added 
subroutines  in  this  program  are  subroutine  MATRIXMULT ,  which 
has  been  described  before,  and  the  matrix  inversion  routine 
INVERSE . 

Subroutine  INVERSE  conducts  the  matrix  inversion  of  the 
eigenvector  matrix  through  a  Gauss-Jordan  routine  (32,*  78). 

The  inversion  is  conducted  by  rotating  the  input  matrix  about 
its  diagonal,  then  normalizing  the  result.  It  is  interesting 
to  note  that  for  the  k  gains  not  to  achieve  an  infinite  value, 
the  controllability  criteria  for  the  system  must  be  satisfied 
for  all  time  increments.  Subroutine  INVERSE  will  stop  the 
execution  of  the  program  if  this  criteria  is  not  met. 
Eigenproblem  Routine 

The  eigenproblem  routine  (Appendix  F)  was  a  converted 
routine  originally  written  in  FORTRAN  4  (17).  This  program 
provides  real  or  complex  eigenvalues  and  eigenvectors  for  a  real 
matrix  in  general  form.  The  program  first  reduces  the  matrix 
to  an  upper  Hessenberg  form  by  means  of  a  Householder  trans¬ 
formation  technique.  Then,  a  unitary  triangular  transformation 
(known  as  a  QR  transformation)  is  conducted  on  the  converted 
Hessenberg  matrix  until  all  elements  of  the  subdiagonal 
converge  to  a  value  which  can  be  approximated  to  equal  zero. 

The  eigenvalues  are  then  determined  from  this  reduced  form. 
Inverse  iteration  is  performed  on  the  upper  Hessenberg  matrix 
until  the  absolute  value  of  the  largest  component  of  the  right 
hand  side  vector  is  greater  than  the  bound  of  100  (a  bound 
established  by  the  numerical  accuracy  performance  of  the 


computer  used).  This  computation  yields  the  eigenvectors  of 
the  matrix.  A  brief  description  of  the  function  of  each 
routine  follows. 

The  executive  routine  of  the  program  performs  the  data 
management  task  of  the  program.  It  first  accepts  the  input 
of  the  monodromy  matrix  as  provided  from  the  program  of 
Appendix  E.  Next,  it  normalizes  the  matrix  through  a  call 
to  subroutine  SCALEE.  Following  this,  it  begins  to  determine 
the  eigenvalues  of  the  normalized  matrix  by  a  call  to  sub¬ 
routine  HESQR.  From  the  eigenvalues,  the  executive  program 
determines  if  the  eigenvectors  will  be  either  real  or  complex, 
and  calls  the  appropriate  subroutine.  Finally,  it  recon¬ 
structs  the  original  matrix  and  converts  the  eigenvectors 
to  this  matrix.  After  the  eigenvectors  are  normalized,  the 
eigenvalues  and  eigenvectors  are  printed  out. 

Subroutine  SCALEE  performs  the  scaling  function  of  the 
original  matrix.  The  array  is  scaled  so  that  the  quotient 
of  the  absolute  sum  of  the  off  diagonal  terms  of  each  column 
and  the  absolute  sum  of  the  off  diagonal  elements  of  each  row 
lies  within  the  bounds  computed  in  the  program.  After  the 
matrix  is  scaled,  it  is  normalized  so  that  the  value  of  the 
Euclidian  norm  is  equal  to  one.  The  scaling  factor  computed 
and  used  in  this  routine  is  stored  for  later  use  in  re¬ 
constructing  the  original  matrix. 

Subroutine  HESQR  determines  the  eigenvalues  of  the  real 
general  matrix.  The  original  matrix  is  converted  to  an  upper 
Hessenberg  form  using  the  Householder  transformation  method. 


The  transformation  process  to  do  this  conversion  uses  a  QR 
approach.  The  end  result  of  this  routine  is  the  indication 
of  whether  the  eigenvalues  were  determined,  the  eigenvalues, 
and  a  flag  to  signify  whether  the  eigenvalues  are  real  or 
complex . 

Subroutine  REALVE  determines  the  real  eigenvectors  of  the 
real  eigenvalues  of  the  matrix.  This  is  performed  using  an 
inverse  iteration  method  which  employs  a  Gaussian  elimination 
technique.  Afterwards,  backsubstitution  is  performed  until 
the  eigenvector  matrix  is  obtained. 

Subroutine  COMPVE  determines  the  complex  eigenvectors  of 
the  complex  eigenvalues  of  the  matrix.  This  routine  splits 
the  complex  eigenvalue  matrix  into  two  matrices  containing  the 
real  parts  and  imaginary  parts  of  the  eigenvalues.  It  then 
uses  an  inverse  iteration  technique,  also  employing  a  Gaussian 
elimination  procedure,  until  the  complex  eigenvectors  are 
determined . 

Subroutine  INPUTT  was  added  to  this  routine  to  handle  the 
input  of  the  monodromy  matrix  to  the  program. 

Feedback  Routine 

The  last  program  used  in  this  analysis  can  be  found  in 
Appendix  J  and  represents  the  system  with  feedback.  The 
structure  of  this  routine  is  very  similar  to  the’  program  of 
Appendix  E  with  the  addition  of  the  feedback  loop  incorporated 
in  subroutine  FCNN.  This  addition  is  found  in  lines  1490-1560 
of  the  program  listing.  The  predetermined  k  gain  value  and 
the  B  matrix  is  inputted  in  lines  1450-1470. 


Step  Size  Determination 

As  the  case  for  all  digital  simulations  of  dynamic  systems, 
the  step  size  selected  for  the  integration  routines  will  have 
a  major  impact  on  the  results  and  the  inferred  stability 
determination  derived  from  these  results .  Because  of  the 
property  of  the  one  sample  time  delay  characteristic  in  all 
digital  simulations,  a  step  size  chosen  too  large  would  infer 
stability  when  in  fact,  stability  does  not  exist  in  the  con¬ 
tinuous  dynamic  system. 

Another  important  consideration  in  any  digital  simulation, 
which  has  a  major  impact  for  the  helicopter  blade  model,  is 
selecting  a  sample  size  so  as  to  capture  all  the  environmental 
factors  which  would  impact  the  stability  of  the  system.  For 
the  helicopter  problem,  the  effects  of  the  various  flow  regions 
as  described  in  Section  II  are  of  primary  concern.  The  width 
of  the  mixed  flow  region  for  the  rotor  is  2.8  degrees  for  the 
conditions  run  for  this  study.  A  step  size  too  large  would 
completely  miss  this  region  and  would  thus  lead  to  erroneous 
results.  For  this  region,  it  is  highly  desirable  to  have 
several  "hits"  in  this  region  to  realize  the  effects  of  mixed 
flow  on  the  flapping  motion  of  the  blade . 

Keeping  these  two  concerns  in  mind,  the  monodromy  matrix 
determination  program  of  Appendix  E  was  run  for  the  various 
azimuth  degree  increments  of  45,30,20,10,5,2,1,  and  0.5 
degrees.  The  resulting  matrices  from  these  runs  were  then 
inputted  into  the  eigenprogram  routine  of  Appendix  F.  The 
resulting  eigenvalues  were  compared  to  determine  convergency 


of  the  roots.  The  results  of  this  analysis  can  be  found  in 
Table  IV . 

As  seen  in  Table  IV,  the  eigenvalues  converge  to  a  constant 
value  at  a  degree  increment  of  one  degree.  Thus,  an  integration 
step  size  of  one  degree  was  selected  for  both  the  monodromy 
matrix  determination  and  the  eigenvector  determination  programs. 
This  selection  led  to  a  program  execution  time  of  41  minutes. 

Because  the  Simpson  integration  routine  used  in  the  Fourier 
series  analysis  programs  uses  data  from  a  one  degree  step  size 
Runga-Kutta  routine,  ten  degree  increment  step  sizes  were  felt 
appropriate  for  these  routines. 

Analysis  Procedure 

The  sequential  use  of  each  program  used  for  this  research 
is  as  follows: 

First,  the  monodromy  matrix  was  determined  using  the  program 
of  Appendix  E,  with  the  appropriate  input  values. 

Next,  the  monodromy  matrix  was  inputted  into  the  eigen- 
problem  program  of  Appendix  F.  This  gave  the  eigenvalues  and 
eigenvectors  of  the  monodromy  matrix.  The  eigenvalues  were 
then  converted  to  the  Poincare'  exponents  to  be  used  in  the 
next  program. 

With  the  information  of  the  Poincare'  exponents,  the 
eigenvectors  of  the  monodromy  matrix  along  with  the  operating 
parameters  used  in  the  first  program  were  used  in  the  program 
of  Appendix  G  to  determine  the  eigenvector  components  for 
every  ten  degrees  of  azimuth  angle. 


Finally,  the  eigenvector  information  was  fed  into  the 
two  Fourierseries  analysis  programs  of  Appendix  H  and  I  to 


determine  the  Fourier  series  representation  of  the  eigenvector 
matrix  and  mode  controllability  matrix. 

Using  the  information  from  all  the  programs  run,  the 
desired  gain  value  for  k  was  determined  for  eigenvalue 
placement.  This  value  of  k  along  with  the  appropriate  B 
matrix  was  inputted  into  the  feedback  control  simulation  of 
Appendix  J  to  verify  this  placement. 
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24. 92593E-02 

63. 20230E-02 

12. 24733E-01 

ONGLE= 

360 

-2 1 . 78727E-02 

43.27003E-02 

95. 39180E-02 

90. 15412E-02 

72 


INCREMENT  OF  RUN  EQUALS  360  STEPS 
FOR  THIS  RUN; 

B=  .$>7  MU=  2.46  GAMMA=  2.00  PSQUARE=  1.00 
POINCARE  EXPONENTS 

OMEGA 1 =-42. 54970E-02  0MEGA2=-62. 33300E-03 

SYSTEM  OF  EIGENUECTORS 

16.5818E-02  31.1379E-02 

-98.61 56E-02  95 . 0286E-02 

A  0  COEFFICIENTS 

-19.6762E-02  -44.8445E-02 

83.71 76E-03  27 . 9506E-03 

B  0  COEFFICIENTS 

. 0000E+00  . 0000E+00 

. 0000E+00  . 0000E+00 

A  1  COEFFICIENTS 

16. 1 161E-02  43.3081E-02 

-59. 7067E-02  52. 9658E-82 

B  1  COEFFICIENTS 

-52. 8504E-02  55.6940E-02 

63 . 7 1 76E-03  -47 . 2890E-02 

A  2  COEFFICIENTS 

58.6173E-03  10.4478E-02 

-42. 6056E-02  26.9492E-02 

B  2  COEFFICIENTS 

-20 . 0555E-02  1 3 . 8087E-02 

-3 1 . 89S9E-03  -2 1 . 7567E-02 

A  3  COEFFICIENTS 

37 . 0986E-03  97 . 6622E-04 

-52. 4787E-03  13. 1302E -02 

B  3  COEFFICIENTS 

43.9555E-03 
-32.8435E-03 
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-12.2395E-03 

-10.6031E-02 


ft  4  COEFFICIENTS 


54 . 9508E-04  -12. 3290E-03 

23 . 7879E-03  -61. 9793E-04 

B  4  COEFFICIENTS 

65.3177E-04  -17.3308E-04 

-24.7809E-03  49.3626E-03 

ft  5  COEFFICIENTS 

1 3 . 8490E-04  -35 . 0863E-04 

1 9 . 5926E-03  11.1 942E-03 

B  5  COEFFICIENTS 

40. 2743E-04  21.8504E-04 

-86 . 2 1 34E-04  1 7 . 4278E-03 


ft  6  COEFFICIENTS 

14.8762E-05  -66.2052E-05 

23. 051 2E-04  -39. 5225E-04 

B  6  COEFFICIENTS 


39 . 1 465E-05  -65 .311 7E-05 

-11. 6032E-04  38 . 0558E-04 

ft  7  COEFFICIENTS 

88.5913E-06  -35.1997E-05 

23 . 4672E-04  36 . 7084E-04 

B  7  COEFFICIENTS 

33 . 7539E-05  53 . 2866E-05 

-74. 4634E-05  25 . 6373E-04 

ft  8  COEFFICIENTS 

1 0 . 9356E-05  83 . 3360E-06 

-26 . 9385E-05  -33 . 3286E-05 

B  8  COEFFICIENTS 

-38 . 7572E-06  -28 . 7435E-06 

-11. 7023E-04  - 1 1 . 4574E-04 
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■v  /.  , 


0  9  COEFFICIENTS 


72 . 1 498E-07  -59 . 5302E-06 

24. 4448E-05  74.3820E-05 


B  9  COEFFICIENTS 


55. 4548E-06  16.5221E-05 

-15. 6874E-05  1 0 . 5438E-04 

010  COEFFICIENTS 


1 0 . 9672E-06  46 . 1 894E-06 

-58. 2987E-06  25. 7 S81E-05 


B10  COEFFICIENTS 

-25. 6176E-06  54.5014E-0? 

-84 . 9673E-05  -11. 0347E-04 

fill  COEFFICIENTS 


-32 . 6572E-06  97 . 0767E-06 

-82. 2738E-05  -11. 0370E-04 

611  COEFFICIENTS 

10.831 8E-05  1 8 . 7380E-05 

-2 1 . 4283E-05  1 0 . 6874E-04 

012  COEFFICIENTS 


-34 . 1 523E-06  25 . 2367E-05 

-80.8118E-05  14. 5163E-04 

B12  COEFFICIENTS 

1 2 . 8032E-05  -20 . 5924E-05 

-56. 1285E-05  90.4627E-05 

013  COEFFICIENTS 

-46 . 406 1 E-05  1 1 . 6634E-04 

-66 . 1 746E-04  -38 . 3395E-04 

B13  COEFFICIENTS 

71 . 7985E-05 
58. 4432E-04 


7  ^ 


13.3460E-04 

-28.3837E-04 


ft  14  COEFFICIENTS 


18.2899E-04 
79. 2725E-84 


41 . 1753E-04 
21.33S7E-04 


B14  COEFFICIENTS 


21.7785E-04 
82. 9825E-04 


-57. 1831E-05 
16. 3374E-03 


ft 15  COEFFICIENTS 


12.3667E-03 
17. 4421E-03 


-32.5420E-04 

-43.8616E-03 


B15  COEFFICIENTS 


40. 8463E-04 
■35. 3442E-03 


14.6429E-03 
-10. 696 IE-03 


ftlt  COEFFICIENTS 


19. 5407E-03 
14. 2029E-02 


-34. 8281E-03 
-89. 8134E-03 


B16  COEFFICIENTS 


■66. 8509E-03 
10.6292E-03 


46. 0046E-03 
-72.5378E-03 


ft 17  COEFFICIENTS 


•53.7198E-03 
19. 9005E-02 


-14.6025E-02 

-17.6603E-92 


B17  COEFFICIENTS 


•17.6170E-02 

21.2431E-03 


18. 5644E-02 
-15. 7610E-02 


ft 18  COEFFICIENTS 


65. 5847E-03 
-27.8947E-03 


14. 9477E-02 
-93. 1 753E-04 


B18  COEFFICIENTS 


-39. 1895E-11 
38.5025E-11 


16.6863E-1 1 
-80.9785E-11 


INCREMENT  OF  RUN  EQUALS  36©  STEPS 


FOR  THIS  RUN;  _  ,  .. 

, 97  MU=  2.40  6AMMA=  4.00  PSQUARE=  1.00 

POINCARE  EXPONENTS 


OMEGA 1 =- 1 0 . 09043E-0 1  0ME6A2=  33.37400E-03 


SVSTEM  OF  EIGENVECTORS 

19. 5071E-02  27.369SE-02 

98.0789E-02  96. 1816E-02 

A  0  COEFFICIENTS 

■30.0984E-02  -52.5723E-02 

30.3711 E-02  -17. 5398E-03 

B  0  COEFFICIENTS 

. 0000E+00  • 0000E+00 

.  0000E+00  ..0000E+00 

A  1  COEFFICIENTS 

92 . 8545E-03  30 . 0678E-02 

-46 . 6463E-02  49. 0Z67E-Q2 

B  1  COEFFICIENTS 

-37 . 2792E-02  48 . 0288E-02 

28 , 3306E-82  -28. 4656E-02 

A  2  COEFFICIENTS 

1 4 . 2285E-02  1 9 . 9223E-02 

-61 . 3334E-02  26.4648E-02 

B  2  COEFFICIENTS 


-23 . 4882E-02  1 2 . 9002E-02 

-47.5636E-03  -39.4158E-02 

A  3  COEFFICIENTS 

87. 9682E-03  57.0997E-03 

-15. 7750E-02  1 4 . 5492E-02 

B  3  COEFFICIENTS 


-23. 0123E-03 
-24.0689E-02 


47. 8377E-03 
-16.9718E-02 


ft  9  COEFFICIENTS 


•10.  7635E-06 
44.0501E-05 


-90. 01 10E- 
12. 4261 E- 


B  9  COEFFICIENTS 


95.7118E-06 

93.6313E-06 


27. 7400E- 
16.2429E- 


ft  10  COEFFICIENTS 


42.5617E-06 

-67.3533E-06 


95. 3395E- 
43. 4979E- 


■06 

■04 


■05 

•04 


■06 

■05 


B10  COEFFICIENTS 


28.6337E-06 

-13.8270E-04 


39. 2950E- 
-14. 9439E' 


fill  COEFFICIENTS 


12.4900E-05 
-20. 5656E-04 


26. 0796E- 
-18.7068E- 


Bll  COEFFICIENTS 


25. 4520E-05 
63. 7637E-05 


33.2638E- 

22.8709E- 


•06 

■04 


■05 

■04 


■05 

■04 


ft 12  COEFFICIENTS 


16. 1 772E-05 
-44.3213E-04 


78.2937E- 
10. 9768E- 


B12  COEFFICIENTS 


72.6581E-05 

-38.4941E-05 


50.5044E- 

36.9516E- 


013  COEFFICIENTS 


-54. 7486E-05 
-17. 7424E-03 


23. 2860E- 
-12. 1545E- 


■05 

■05 


■06 

-04 


■04 

■03 


B13  COEFFICIENTS 


36.0658E-64 

-63.8731E-04 


23.9805E- 
1 1 . 7349E- 


■04 

■03 


Hl4  coefficients 


INCREMENT  OF  RUN  EQUALS  360  STEPS 


»  L 


FOR  THIS  RUN; 

B=  .97  HU*  2.40 


PSQUARE*  1.00 


6AMMA=  6.00 
POINCARE  EXPONENTS 

OMEGA 1 =- 1 5 . 49838E-0 1  0ME6A2-  86.369O0E-O3 

SVSTEM  OF  EIGENVECTORS 


21 . 2337E-02 
-97.7197E-02 


35. 1497E-02 
93.6189E-02 


u  Pi  COEFFICIENTS 


-32. 2339E-02 
49 . 9578E-02 


-54.3168E-02 

-46.8978E-03 


B  0  COEFFICIENTS 

. 00O0E+00  ■ OOOOE+OO 

. 0000E+0Q  . C0OOE+00 


Pl  1  COEFFICIENTS 


21 .6823E-03 
-36. 2954E-02 


26.6427E-02 


49. 9360E-O2 
fi  I  COEFFICIENTS 


-32 . 9368E-02 
48. 877SE-02 

A  2 


47.6285E-02 

-22.5302E-02 


COEFFICIENTS 


17.6584E-02 

-75.3277E-02 


26.9833E-02 

22.0735E-O2 


B  2  COEFFICIENTS 


-23. 980OE-O2 
18.4834E-03 


98.7114E-03 
-53. 1 168E-02 


A  3  COEFFICIENTS 


13.3284E-02 

-28.7070E-02 


98. 71 15E-03 
89. 0562E-03 


B  3  COEFFICIENTS 


-26.8498E-03 

-35.8259E-02 


26. 8O68E-03 
-29. 3855E-02 
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AD-A154  468  APPLICATION  OF  FLOQUET  THEORV  TO  HELICOPTER  BLADE 
FLAPPING  STABILITV(U)  AIR  FORCE  INST  OF  TECH 
HRIGHT-PATTERSON  AFB  OH  SCHOOL  OF  ENGINEERING 
UNCLASSIFIED  J  K  HARCH  DEC  84  AFIT/GAE/AA/84D-12  F/G  1/2 


Pi  4  COEFFICIENTS 


42 . 1 554E-03  -11. 8375E-04 

33 . 9278E-03  37 . 9430E-03 


B  4  COEFFICIENTS 

24 . 8203E-03  95 . 1 989E-04 

-20 . 7 1 33E-02  53 . 9532E-04 

ft  5  COEFFICIENTS 

3 1 . 9844E-04  -65 . 2583E-04 

93 . 4039E-03  74 . 7848E-03 

B  5  COEFFICIENTS 

19.6558E-03  15.0479E-03 

-46.4943E-03  33.8744E-03 

ft  6  COEFFICIENTS 

-19. 8555E-04  -29 . 7538E-04 

35.5153E-03  17.9749E-03 

B  6  COEFFICIENTS 

54 . 2028E-04  30. 6696E-04 

33. 0717E-04  1 7 . 5934E-03 

PI  7  COEFFICIENTS 

-14.2159E-04  -14. 1U4E-04 

11.5912E-03  13.9712E-03 

B  7  COEFFICIENTS 


13.3851E-04  20.5031E-04 

78.3281E-04  10.0977E-03 

ft  8  COEFFICIENTS 

-13.0515E-05  -33.5245E-06 

12.7509E-04  14.7852E-04 


B  8  COEFFICIENTS 

1 6 . 3824E-65  24 . 6929E-05 

22 . 2372E-05  -82 . 8043E-05 
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•A 


A  9  COEFFICIENTS 


-46. 5S16E-06 
46. 4886E-05 


-11.0951E-05 

19.8482E-04 


B  9  COEFFICIENTS 


88.2182E-06 
70. 39-54E-05 


44.4661E-05 

20.3948E-04 


B10  COEFFICIENTS 


13.0906E-05 

-20.0284E-05 


18.3499E-05 

25.8491E-05 


B10  COEFFICIENTS 


90.3189E-06 
-83. 1601E-05 


15.5450E-05 
-19. 465 2E “04 


All  COEFFICIENTS 


48.0838E-05 
-39. 1808E-04 


46.441 6E-05 
-42.6440E-84 


81 1  COEFFICIENTS 


44. 1368E-05 
25.4166E-04 


71.6457E-05 

34.3360^-04 


A12  COEFFICIENTS 


69.3890E-05 
-11. 7335E-03 


10. 7174E-04 
-56.8159E-04 


812  COEFFICIENTS 


18. 1969E-04 
72.4322E-05 


10.4698E-04 

49.4734E-04 


A13  COEFFICIENTS 


-10.5981E-04 
-31 . 3044E-03 


21.8299E-04 
-25. 1487E-03 


B13  COEFFICIENTS 


65.3715E-04 

-15.5567E-03 


49. 9286E-04 
11. 1807E-03 


014 

COEFFICIENTS 

14.0456E-03 

11.2766E-03 

41.5557E-05 
-12. 5759E-03 

614 

COEFFICIENTS 

82.7694E-04 
■69. 1282E-03 

31.7815E-04 
14. 8758E-04 

015 

COEFFICIENTS 

•44. 4240E-03 
95.5889E-03 

-32. 8994E-03 
-29.9151E-03 

BIS 

COEFFICIENTS 

■89.5944E-04 
■11. 9456E-02 

89. 1404E-04 
-98. 0347E-03 

016 

COEFFICIENTS 

58.8642E-03 
25. 1095E-02 

-89. 9490E-03 
-73.5972E-03 

616 

COEFFICIENTS 

79. 9327E-03 
61 . 6862E-04 

32.9036E-63 

-17.7106E-02 

017 

COEFFICIENTS 

■72.2501E-04 

12.0953E-02 

-88.8079E-03 

-16.6585E-02 

817 

COEFFICIENTS 

■10.9792E-02 

16.2914E-02 

15.8754E-02 
-75. 1278E-03 

018 

COEFFICIENTS 

10.7442E-02 

16.6533E-02 

18. 1045E-02 
15.5927E-03 

618 

COEFFICIENTS 

■48.8196E-11 
94.8963E-1 1 

-91.2150E-12 
-90.8258E-1 1 

INCREMENT  OF  RUN  EQUALS  36  STEPS 


FOR  THIS  RUN; 

B=  .97  MU=  2.40  GAMMA=  8.00  P 
POINCARE  EXPONENTS 

OMEGA 1 =-20. 36572E-01  0ME6A2=  95 

SVSTEM  OF  EIGENVECTORS 

20.5712E-02  43.25O9E-02 

•97.8613E-02  90. 1630E-02 

A  O  COEFFICIENTS 

11.8772E-01  10. 2515E-01 

86 . 2073E-02  1 3 . 6359E-0 1 

B  0  COEFFICIENTS 

. 0000E+00  . 0000E+00 

. 0000E+00  • 0000E+00 

A  1  COEFFICIENTS 

-21.5499E-02  46.5291E-02 

1 9 . 6603E-02  11. 3453E-01 

B  1  COEFFICIENTS 

-28.9105E-02  13.4273E-02 

86. 4842E-03  -21 . 5431E-02 

A  2  COEFFICIENTS 

1 4 . 8039E-02  84 . 5943E-02 

95.8541E-02  70.8016E-02 

B  2  COEFFICIENTS 

42 . 7057E-03  36 . 8620E-03 

39 . 856 1 E-02  -23 . 6625E-02 

A  3  COEFFICIENTS 

23. 1093E-02  68. 2755E-02 

57. 5519E-02  88. 9007E-02 

B  3  COEFFICIENTS 


78.6651E-04 
-20. 2031 E-02 


-89. 6951E-03 
-58.3766E-02 


Pi  4  COEFFICIENTS 


81.0408E-03 

59.7315E-03 


42. 7697E-02 
-13.5355E-02 


B  4  COEFFICIENTS 


-22. 0729E-02 
-40.0155E-02 


-29. 1448E-02 
-17. 2588E-02 


ft  5  COEFFICIENTS 


10.8282E-02 

10.8934E-02 


39. 9897E-02 
37 . 3628E-02 


B  5  COEFFICIENTS 


84. 7035E-03 
75.8592E-02 


15. 2980E-02 
20. 3044E-02 


ft  6  COEFFICIENTS 


27.9636E-02 

73.8246E-02 


49. 7259E-02 
50. 151 IE-02 


B  6  COEFFICIENTS 


-74.8492E-02 

-33.8013E-04 


-67.5640E-02 

-15.4892E-02 


0  7  COEFFICIENTS 


-68.7952E-02 

-48.5174E-03 


-56.6376E-02 
-26. 1896E-02 


B  7  COEFFICIENTS 


41.0381E-02 
65. 7802E-02 


57 . 0929E-02 
44.0293E-02 


ft  8  COEFFICIENTS 

7 1 . 6266E-03  l 3 . 5705E-02 

4 1 . 6068E-02  1 7 . 5857E-02 

B  8  COEFFICIENTS 


13. 9513E-02 
15.3235E-02 


35.3739E-02 

56.4575E-02 


U  9  COEFFICIENTS 


45.0195E-03 
-48 . 6608E-03 


S5.3710E-03 

-12.3636E-03 


COEFFICIENTS 


84.36S9E-04 

47.7471E-02 


29.6658E-02 
17. 1040E-02 


COEFFICIENTS 


-48.3488E-02 
-31 . 0404E-02 


-51. 1885E-02 
-35. 9755E-02 


COEFFICIENTS 


23.2088E-04 

71.4163E-02 


31.3986E-02 
92. 1507E-02 


COEFFICIENTS 


23.0722E-03 

88.4364E-02 


21.5094E-03 

37.7143E-02 


COEFFICIENTS 


34.8169E-02 

72.0456E-02 


72.0564E-02 

59.0006E-02 


COEFFICIENTS 


-35. 242OE-02 
12.2056E-02 


-38.6831E-02 
77 . 5380E-03 


COEFFICIENTS 


-26.8505E-03 
-19. 1759E-02 


38. 7129E-02 
34.9878E-02 


COEFFICIENTS 


12.8885E-02 
-11. 9757E-02 


80. 7772E-03 
-97.091 IE-03 


COEFFICIENTS 


10. 1510E-01 
92.2498E-02 


14.9400E-01 

10.7946E-01 


ftl4  COEFFICIENTS 


20. 6526E-02 
13. 7431E-82 


16.5659E-02 

33.9748E-02 


B14  COEFFICIENTS 


81.9684E-03 

10.6607E-02 


30. 7876E-02 
43.64S2E-02 


015  COEFFICIENTS 


57.4714E-03 
45. 0695E-03 


-11. 5097E-02 
-17.3672E-02 


B15  COEFFICIENTS 


78.3321E-03 

45.4060E-02 


46.9790E-02 

22.6939E-02 


ft 16  COEFFICIENTS 


-26. 1592E-02 
-32.6501E-02 


-42. 3541 £-02 
-55.6018E-03 


B16  COEFFICIENTS 


54. 4707E-02 
47.68S5E-02 


94.7728E-02 

58.0904E-02 


017  COEFFICIENTS 


63.2308E-02 

16.9784E-03 


51 . S229E-02 
-13. 0551E-02 


B17  COEFFICIENTS 


36.0755E-03 

-41.7061E-02 


49.6161E-02 

-19.7202E-02 


ftlS  COEFFICIENTS 


-77.5633E-03 
-99. 7068E-02 


74.5561E-03 

-45.6649E-02 


B18  COEFFICIENTS 


54.3097E-03 

20.8531E-02 


32. 7916E-02 
20. 1257E-02 


INCREMENT  OF  RUN  EQUALS  360  STEPS 


FOR  THIS  RUN; 

B=  .97  MU=  2.40 


PSQUARE=  1.00 


SAMMA=10. 00 
POINCARE  EXPONENTS 

OMEGA 1 =-25. 12578E-01  0MEGA2=  81.90400E-03 

SYSTEM  OF  EIGENUECTORS 


-18.3184E-02 

98.3079E-02 


50. 7958E-02 
86. 1382E-02 


A  0  COEFFICIENTS 


33.9228E-02 

-83.4099E-02 


-56. 6455E-02 
-46. 3378E-03 


8  0  COEFFICIENTS 


. 080OE+00  • 0000E+0® 

. 0800E+00  • 0000E+00 


A  1  COEFFICIENTS 

94.0188E-03  24.6571E-02 

89. 0624E-03  54.6691E-02 

B  1  COEFFICIENTS 

30 . 7857E-02  52 . 637 1 E-02 

-86 . 5526E-02  -20.3481 E-02 

A  2  COEFFICIENTS 


-15.6284E-02 

93.9369E-02 


37.4191E-02 
15. 6794E-02 


8  2  COEFFICIENTS 


26. 4231 E-02 
-35. 1361 E-02 


63. 0430E-03 
-74.3271E-02 


A  3  COEFFICIENTS 


-18.2564E-02 

63.5598E-02 


15. 2363E-02 
-83.6715E-03 


8  3  COEFFICIENTS 


52.8993E-03 
41 .4744E-02 


-32. 1217E-03 
-45.9820E-02 
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ft  4  COEFFICIENTS 


86. 1619E-03  16. 1701E-03 

70. 9883E-03  16.2254E-05 


6  4  COEFFICIENTS 


40 . 9324E-03  -29 . 5246E-05 

44. 7398E-02  -65.0169E-03 

ft  5  COEFFICIENTS 

•14.3148E-03  19.1 124E-04 

•15.7159E-02  13. 3819E-02 

B  5  COEFFICIENTS 

-42 . 2562E-03  26 . 6900E-03 

1 7 . 7584E-02  -76 . 3040E-04 

fl  6  COEFFICIENTS 


56.111 2E-04  43 . 8628E-05 

-10. 8998E-02  65 . 982 1 E-03 

B  6  COEFFICIENTS 

-18. 8645E-03  1 1 . 0360E-03 

13. 4928E-03  -27.2000E-04 

ft  7  COEFFICIENTS 

58 . 0747E-04  -10. 2628E-04 

-40.21 35E-03  4 1 . 8699E-03 

B  7  COEFFICIENTS 

-62 . 4839E-04  60 . 7098E-64 

-25.51 05E-03  73 . 7344E-04 

ft  8  COEFFICIENTS 


24.4285E-04  -20.3257E-06 

-75.5173E-04  11.7410E-03 

B  8  COEFFICIENTS 

-24. 4200E-04  16.5948E-04 

-14.41 11E-03  -15. 5089E-04 


ft  9  COEFFICIENTS 


88. 2273E-05 
•34.5439E-05 


-15.8638E-05 

51.8213E-64 


B  9  COEFFICIENTS 


•17. 4481E-04 
-77.7139E-04 


1 1 . 5727E-04 
29.9064E-04 


ft 1 0  COEFFICIENTS 


-14.9817E-05 
22. 1497E-04 


29.3766E-05 

-16.2925E-04 


B10  COEFFICIENTS 


-18. 8793E-04 
-59.4440E-04 


77. 8146E-05 
-32. 7850E-04 


fill  COEFFICIENTS 


-13.6594E-04 
13. 1050E-03 


39.0955E-05 
-13. 1224E-03 


Bll  COEFFICIENTS 


-29.9311E-04 

-93.0120E-04 


20.9393E-04 
19. 5546E-04 


ftl2  COEFFICIENTS 


-13. 7497E-04 
35. 7185E-03 


52. 9246E-07 
-21 . 5120E-03 


B12  COEFFICIENTS 


-  70. 7485E-04 
40.8066E-04 


37. 1427E-04 
-26. 4279E-04 


ftl3  COEFFICIENTS 


52.3338E-04 

52.0155E-03 


-59.9131E-05 

-45.0293E-03 


B13  COEFFICIENTS 


-14. 7054E-03 
58. 7443E-03 


88.5219E-04 

-30.3782E-04 


AH  COEFFICIENTS 


2'?.  1538E-03 
24. 1669E-03 


-53. 4949E-04 
-86. 0706E-06 


B14  COEFFICIENTS 


14. 1407E-03 
14.8762E-02 


-10. 5155E-05 
-22. 2716E-03 


A15  COEFFICIENTS 


61 . 2629E-03 
21. 2259E-02 


-50. 7777E-03 
27. 4485E-03 


B15  COEFFICIENTS 


17. 2748E-03 
13. 7989E-02 


-10. 7488E-03 
-15. 3485E-02 


PI 6  COEFFICIENTS 


52. 4931E-03 
31 . 3567E-02 


-12.4739E-02 
-52. 3 904E-03 


B16  COEFFICIENTS 


87. 8363E-03 
11. 7333E-02 


21.0055E-03 
-24. 7853E-02 


A17  COEFFICIENTS 


30.9550E-03 
30. 1174E-03 


-82. 1938E-03 
-18.24 83E-02 


B17  COEFFICIENTS 


10.2235E-02 

-23.3593E-02 


17. 5442E-02 
-67. S751E-03 


A18  COEFFICIENTS 


-11.2639E-02 
27. 7605E-02 


13. 3797E-02 
15. 3229E-03 


B18  COEFFICIENTS 


47. 2773E-1 1 
-17. 2339E-10 


-25.5162E-1 1 
-89.0636E-11 


If 


BB 


r 


p  i 


■  « 


r.- 

. 


»  i 


INCREMENT  OF  RUN  EQUALS  360  STEP'i 


FOR  THIS  RUN; 

B=  .97  MU®  2.40 


PSQUARE®  1.00 


GAMMA® 12. 00 
POINCARE  EXPONENTS 

OMEGA 1 =-29 . 86568E-0 1  0MEGA2=  56. 85900E-0O 

SVSTEM  OF  EIGENUECTORS 


-15. 1443E-02 
98. 8466E-02 


57. 8038E-02 
81 . 6010E-02 


A  0  COEFFICIENTS 


15.0218E-02 

-31.2415E-01 


-59. 1439E-02 
-33. 551 6E -03 


G  0  COEFFICIENTS 


. 0000E+00 
. 0000E+00 


. 8000E+08 
. 0000E+00 


A  1  COEFFICIENTS 


-16.6921E-02 
-19. 8158E-01 


24 . 4537E-G2 
58 . 0557E-02 


B  1  COEFFICIENTS 


19.4244E-02 
-41. 1836E-02 


56. 6499E-02 
-21.2354E-02 


A  2  COEFFICIENTS 


-59.8723E-02 
-42. 7799E-02 


42. 3082E-02 
13. 472SE-02 


B  2  COEFFICIENTS 


22.7456E-82 

52.2407E-02 


55.2921E-03 
-84. 3089E—02 


A  3  COEFFICIENTS 


-70. 4334E-02 
-48. 3020E-03 


17. 2261 £-02 
-17. 7768E-02 


B  3  COEFFICIENTS 


17.3728E-02 
16. 0456E-01 


-62. 6099E-03 
-52. 0484E-02 
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fmm.  V 


Pi  4 

COEFFICIENTS 

54. 9770E-02 
29.6615E-02 

21.2675E-03 

-36.9130E-03 

B  4 

COEFFICIENTS 

17.9551E-02 

16.8450E-01 

-95. 4169E-04 
-86.0159E-03 

PI  5 

COEFFICIENTS 

36.9613E-02 

43.9165E-02 

70. 0I40E-04 
14. 9868E-02 

B  5 

COEFFICIENTS 

22.0589E-02 

12.2973E-01 

29. 8431E-03 
-33. 7015E-03 

fi  6 

COEFFICIENTS 

25.8603E-02 
32. 7618E-02 

40. 7577E-04 
86. 8063E-03 

B  6 

COEFFICIENTS 

25.5376E-02 

85.8603E-02 

14. 4908E-03 
-24.8821E-03 

ft  7 

COEFFICIENTS 

19.5167E-02 

17.5949E-02 

27. 9879E-05 
57. 8227E-03 

B  7 

COEFFICIENTS 

26.5199E-02 
63. 7266E-02 

83. 7907E-04 
-19. 8926E-04 

Pi  8 

COEFFICIENTS 

15.2316E-02 

•73.5266E-03 

36. 1 309E-05 
19.6705E-03 

B  8 

COEFFICIENTS 

25.8109E-02 
62. 3278E-02 

27.7880E-04 

-50.0961E-04 

ft 14  COEFFICIENTS 


91.7290E-03 

13.7644E-02 


-70. 

12. 


3685E-04 

2064E-03 


B14  COEFFICIENTS 


12.7416E-02 

70.2855E-02 


-31. 

-29. 


9444E-04 

4447E-03 


fil5  COEFFICIENTS 


14.8245E-02 

54.6199E-03 


-57. 

58. 


4064E-03 

6905E-03 


B15  COEFFICIENTS 


10.8113E-02 

83.8075E-02 


-20. 

-17. 


9230E-03 

3790E-02 


ftl6  COEFFICIENTS 


1 1 . 6430E-02 
18.4950E-02 


-14 

-45 


1038E-02 

1034E-03 


B16  COEFFICIENTS 


10.9054E-02 
24. 1557E-02 


18 

-28 


.4163E-03 
.  1 156E-02 


fll7  COEFFICIENTS 


-25.5512E-03 

70.3416E-02 


-81 

-19 


.51S5E-03 
,  3843E-02 


B17  COEFFICIENTS 


81.2972E-03 

-10.3954E-02 


18 

-70 


,  8814E-02 
,  8469E-03 


fll8  COEFFICIENTS 


-13.0619E-02 

10.8443E-01 


19 

11 


,  7119E-02 
,  0081E-03 


B18  COEFFICIENTS 


70.5010E-11 

11.2050E-10 


-32 

-91 


,  0229E-11 
,  3379E-11 


RQ 


INCREMENT  OF  RUN  EQUALS  366  STEPS 


FOR  THIS  RUN; 

B~  .9?  MU=  2.40 


PSQUARE*  1.00 


GAMMA* 14. 00 
POINCARE  EXPONENTS 

OMEGA1=-30.31?60E-01  0MEGA2®  25.30000E-03 

SVSTEM  OF  EIGENUECTORS 


-11.3653E-02 
‘  99.3521E-02 


64 . 3938E-02 
76. 5077E-02 


Ti 

O 


,  4 

r 


A  0  COEFFICIENTS 


-28.6256E-02 

-54.5366E-01 


-62.4677E-02 
-15. 6881E-03 


B  0  COEFFICIENTS 


. 0000E+00 
.0000E+00 


. 0000E+00 
.0000E+00 


r-- 

> 


A  1  COEFFICIENTS 


-62.0704E-02 
-44. 7871E-01 


24.4110E-02 
61 . 9486E-02 


B  1  COEFFICIENTS 


-41.6308E-03 

75.1732E-02 


61.3108E-02 
-22. 8636E-02 


A  2  COEFFICIENTS 


-11.6490E-01 

-25.6014E-01 


47. 2723E-02 
1 1 . 6321E-02 


B  2  COEFFICIENTS 


10.9609E-02 
20. 1084E-01 


52. 1135E-03 
-94.4216E-02 


A  3  COEFFICIENTS 


|L- 


-13.3823E-01 

-11.8021E-01 


19.0480E-02 
-27. 3447E-02 


B  3  COEFFICIENTS 


35. 3525E-02 
29. 7245E-01 


-92.8687E-03 

-57.3983E-02 
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*‘.A  «,*  *, 

Sm  ■  •-  U  mJJL* 


A.  t. 


'  V*  •  .  •* .  , 


PI  4  COEFFICIENTS 


11.4031E-01 

68.2079E-02 


24. 1271E-03 
-79.7181E-03 


COEFFICIENTS 


53.5045E-02 

29.9723E-01 


-20. 1034E-03 
-97.5216E-03 


COEFFICIENTS 


85.0163E-02 
57. 2038E-02 


1 1 . 6965E-03 
15.9057E-02 


COEFFICIENTS 


62. 1112E-02 
24.7227E-01 


31.6922E-03 

-58.2183E-03 


COEFFICIENTS 


-62.3884E-02 

-44.8369E-92 


84. 0106E-04 
10. 311 IE-82 


COEFFICIENTS 


64. 7382E-82 
19.5639E-01 


17.2363E-03 
-5 1 . 5006E-03 


COEFFICIENTS 


•47.0747E-02 
■27 . 2968E-02 


22.6032E-04 
72. 4014E-03 


COEFFICIENTS 


64.0689E-02 

16.3363E-01 


10.5232E-03 
-16. 3190E-03 


COEFFICIENTS 


•36.0972E-02 

■11.2331E-02 


10.8893E-04 
28. 1818E-03 


COEFFICIENTS 


61.5902E-02 
14. 5145E-01 


40. 2642E-04 
-11. 5496E-03 


0  9  COEFFICIENTS 


27.2806E-02 
88. 1660E-04 


-18.2606E-0? 
1 1 . 6647E-03 


B  9  COEFFICIENTS 


58.2695E-02 

13.2970E-01 


25. 9673E-04 
18. 7921E-05 


010  COEFFICIENTS 


19.5609E-02 

10.4976E-02 


93. 3743E-06 
-34 . 0530E-04 


B10  COEFFICIENTS 


54 . 6226E-02 
12.3589E-01 


19.3802E-04 
-82. 5379E-04 


Oil  COEFFICIENTS 


•12. 1975E-02 
19.3931E-02 


-63.3059E-05 

-22.0698E-03 


BU  COEFFICIENTS 


56.7416E-02 

11.7277E-01 


36.8685E-04 

-66.8696E-04 


012  COEFFICIENTS 


-42. 3835E-03 
27.4877E-02 


-25.6644E-04 
-33. 4701E-03 


B12  COEFFICIENTS 


46.3925E-02 
1 1.6697E-01 


58. 1242E-04 
-19. 8441E-03 


013  COEFFICIENTS 


54.8264E-03 

32.9964E-02 


-38. 2018E-04 
-53. 6207E-03 


B13  COEFFICIENTS 


41. 1004E-02 
12.3573E-01 


10.4990E-03 

-20.4115E-03 


014  COEFFICIENTS 


16. 7868E-02 
37.4438E-02 


-79. 7684E-04 
26. 4079E-03 


814  COEFFICIENT'; 


33.9602E-02 
13. 1674E-01 


-67 . 2289E-04 
-33. 4724E-03 


015  COEFFICIENTS 


24.5662E-02 

54.4406E-02 


-63. 4748E-03 
90.4510E-03 


B15  COEFFICIENTS 


23.7635E-02 

12.2218E-01 


-3 1 . 0220E-03 
-19. 1721E-02 


016  COEFFICIENTS 


19.5892E-02 

10.0601E-81 


-15. 7587E-02 
-39.0458E-03 


B16  COEFFICIENTS 


11.5841E-02 
82. 1200E-02 


17. 3503E-03 
-31.4901E-02 


017  COEFFICIENTS 


19. 1216E-03 
16. 4605E-01 


-81 . 3795E-03 
-20. 6899E-02 


B17  COEFFICIENTS 


25.6071E-03 

32.5043E-02 


20.4346E-02 

-76.2910E-03 


018  COEFFICIENTS 


-90.8445E-03 

19.7112E-01 


20. 8192E-02 
49.9261E-04 


B18  COEFFICIENTS 


91.6327E-11 
58. 7632E-10 


-38.2674E-11 

-95.2217E-11 


INCREMENT  OF  RUN  EQUALS  360  “STEPS 


FOR  THIS  RUN; 

B=  .97  MU=  2.40  GAMMA=16.00  PSQUARE*  1.00 

POINCARE  EXPONENTS 

OMEGA 1 *-33 . 36009E-0 1  0MEGA2*- 1 0 . 56400E-03  • 

SYSTEM  OF  EIGENUECTORS 

-70 . 7829E-03  70 . 5935E-02 

99. 7492E-02  70 . 8276E-02 

A  0  COEFFICIENTS 

-32.  7740E-01  -66. 3098E-02 

-22.4760E+00  71.4889E-04 

8  0  COEFFICIENTS 

. 0000E+00  . 0000E+00 

. 0000E+00  . 0000E+00 

A  1  COEFFICIENTS 

-4 1 . 6939E-0 1  24 . 3642E-02 

-20 . 0876E+00  66.11 53E-02 

B  1  COEFFICIENTS 

-60.2491 E-02  66 . 3488E-02 

62. 8587E-01  -25. 0699E-02 

A  2  COEFFICIENTS 

-56.4746E-01  52.2731E-02 

-14.3530E+00  10.0159E-02 

B  2  COEFFICIENTS 

77. 6956E-03  52 . 7564E-03 

1 1 . 1085E+00  -10.4623E-01 

A  3  COEFFICIENTS 

-60 . 5523E-0 1  20 . 7920E-02 

-84 . 1 436E-0 1  -36.931 5E-02 

B  3  COEFFICIENTS 

15.5215E-01  -12. 2506E-02 

13. 1475E+00  -62.2716E-0: 


ft  4 

COEFFICIENTS 

-53.0780E-01 

-45.6764E-01 

25 . 2434E-03 
-12. 5413E-02 

B  4 

COEFFICIENTS 

27.3656E-01 

12.4135E+00 

-31 . 3130E-03 
-10. 1257E-02 

ft  5 

COEFFICIENTS 

-41.9517E-01 

-26.0204E-01 

15. 7361E-03 
16. 3169E-02 

B  5 

COEFFICIENTS 

32.8337E-01 
10. 5793E+00 

32.6071E-03 
-79. 7277E-03 

ft  6 

COEFFICIENTS 

-32. 2891E-01 
-16.2615E-01 

13.0426E-03 
1 1 . 4796E-02 

B  6 

COEFFICIENTS 

34. 0446E-01 
89.4790E-01 

19.2846E-03 
-80 . 306 1 E-03 

ft  7 

COEFFICIENTS 

-24. 9460E-01 
-90. 8492E-02 

47. 5139E-04 
84 . 3722E-03 

B  7 

COEFFICIENTS 

33.3579E-01 
78. 1620E-01 

12. 3493E-03 
-34. 5395E-03 

ft  8 

COEFFICIENTS 

-1 9.2706E-01 
-29. 1518E-02 

21 . 4708E-04 
36. 2859E-03 

B  8 

COEFFICIENTS 

31.9267E-01 

70.3433E-01 

52. 8384E-04 
-21.2346E-03 

0  9  COEFFICIENTS 


-14.6052E-01 

22.8217E-02 


24. 1306E-05 
15.7427E-03 


B  9  COEFFICIENTS 


30.2022E-01 

64.4257E-0I 


35. 0538E-04 
-42. 6224E-04 


010  COEFFICIENTS 


-10.5143E-01 

66.6679E-02 


-97.6937E-06 
-33 . 2802E-04 


B10  COEFFICIENTS 


28. 3583E-01 
59. 7441E-01 


26.3691E-04 

-13.0927E-03 


Oil  COEFFICIENTS 


-66. 6088E-02 
10. 5208E-01 


-14. 0386E-04 
-24.9478E-03 


Bll  COEFFICIENTS 


26.4281E-01 

56.3059E-01 


43.9581E-04 

-13.3980E-03 


012  COEFFICIENTS 


-26.9831E-02 

14.0762E-01 


-40.6537E-04 

-36.9449E-03 


B12  COEFFICIENTS 


24. 3018E-01 
54.5844E-01 


65.2893E-04 

-30.0009E-03 


013  COEFFICIENTS 


16.8220E-02 

18.0370E-01 


-51.4264E-04 

-55.0155E-03 


B13  COEFFICIENTS 


21.6091 E-01 
55.0376E-01 


10. 7987E-03 
-27. 8827E-03 
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A14  COEFFICIENTS 


Appendix  D. 

Mode  Controllability  Matrices  Element 
Fourier  Series  Coefficients 
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INCREMENT  OF  RUN  EQUALS  360  STEPS 
FOR  THIS  RUN; 

B=  .9?  MU=  2.40  6AMMA=  2.00  PSQUfiRE®  1 . 00 
POINCARE  EXPONENTS 

OMEGA 1 =-42. 54970E-02  0ME6A2=-62. 33300E-03 

SVSTEM  OF  EIGENUECTORS 

16. 581 8E -02  31 . 1379E-82 

-98 . 6 1 56E-02  95 . 0286E-02 

A  O  COEFFICIENTS 

1 7 . 2936E-0 1  -67.1 089E-02 

B  0  COEFFICIENTS 

. 0000E+O0  . 000OE+O0 

A  1  COEFFICIENTS 

-13. 5909E-01  64. 357 IE-02 

B  1  COEFFICIENTS 

-15. 1 lOOE-Ol  -15. 2794E-01 

A  2  COEFFICIENTS 

-98 . O245E-03  80 . 0865E-03 

B  2  COEFFICIENTS 

-29 . 8593E-02  -35 . 4408E-02 

A  3  COEFFICIENTS 

-12. 5109E-O2  22. 6213E-04 

B  3  COEFFICIENTS 

52 . 5992E-04  - 1 4 . 1 1 99E-03 


H  4  COEFFICIENTS 
44 . 7952E-03  -32 . 3929E-03 

B  4  COEFFICIENTS 
38. 6 335 E -03  28. 9450E— 03 

ft  5  COEFFICIENTS 
-71. 6853E-04  -60 . 49 1 4E-05 

B  5  COEFFICIENTS 
1 8 . 0087E-04  9 1 . 9358E-04 

Pi  6  COEFFICIENTS 
88 . 9999E-04  -29.9810E-04 

B  6  COEFFICIENTS 
58. 921 2E -04  36.3121E-04 

fl  7  COEFFICIENTS 
-15.1 0S6E-04  1 0 . 4907E-04 

B  7  COEFFICIENTS 
-29 . 1 683E-04  23 . 4895E-04 

Pi  8  COEFFICIENTS 
16.61 53E-04  -35 . 3448E-05 

B  8  COEFFICIENTS 
25. 3774E-04  65.4536E-05 

Pi  9  COEFFICIENTS 
1 3 . 5823E-05  35 . 6 1 50E-05 

B  9  COEFFICIENTS 
-20 . 4539E-04  12.31 73E-04 
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010  COEFFICIENTS 


■47.91 26E-05  30 .  7768E-05 

B10  COEFFICIENTS 
1 6 . 2466E-04  28 . 2642E-05 

Fill  COEFFICIENTS 
89.6143E-05  -12.7433E-05 

Bll  COEFFICIENTS 
■15. 833 1 £-04  10.  7334E-04 

012  COEFFICIENTS 
•31. 8666E-04  1 2 . 0797E-04 

B12  COEFFICIENTS 
2 1 . 6224E-04  11.541 0E-04 

013  COEFFICIENTS 
26 . 6464E-04  28 . 4903E-05 

B13  COEFFICIENTS 
4 1 . 0272E-05  31.11 29E-04 

014  COEFFICIENTS 
•15. 1 151E-03  10.9167E-03 

B14  COEFFICIENTS 
1 2 . 8465E-03  95.8140E-04 

015  COEFFICIENTS 
41.8141E-03  -74.6242E-05 

B15  COEFFICIENTS 
17.3734E-04  -47.3791E-04 


ft  16-  COEFFICIENTS 


32.5848E-03  -26.6601E-03 

B16  COEFFICIENTS 
■99.5927E-03  -11.8176E-02 

ftl?  COEFFICIENTS 
45 . 3038E-02  -2 1 . 4549E-02 

B17  COEFFICIENTS 
•50 . 3656E-02  -50. 9333E-02 

ft 18  COEFFICIENTS 
•57 . 6503E-02  22 . 3700E-02 

BIS  COEFFICIENTS 
-21.6881E-11  -88.3967E-1 1 


INCREMENT  OF  RUN  EQUALS  360  STEPS 


B= 


FOR  THIS  RUN:  4  , 

.97  MU=  2.40  SAMMA=  4.00  F‘SQUARE=  1.6 

POINCARE  EXPONENTS 

OMEGA 1=-10.09043E-01  0MEGA2=  33.37400E-03 

SVSTEM  OF  EIGENUECTORS 


19.5071E-02  27 . 3698E-02 

-98. 0789E-02  96. 1816E-02 


A  0  COEFFICIENTS 
30 . 1 820E-0 1  -13. 7756E-0 1 

B  0  COEFFICIENTS 
. 0000E+00  • 0000E+00 

A  1  COEFFICIENTS 
-15.1 648E-0 1  ?8. 6900E-02 

B  1  COEFFICIENTS 
-17. 1629E-01  -14.0303E-01 

A  2  COEFFICIENTS 


-32. 9207E-02  34.7404E-02 

B  2  COEFFICIENTS 
-38 . 3547E-02  -42 . 5539E-02 

A  3  COEFFICIENTS 


-38 . 9096E-02  75 . 54 1 2E-03 

B  3  COEFFICIENTS 
94.6384 E-03  -26.8808E-05 


ft  4  COEFFICIENTS 


10.9579E-02  -80.3465E-03 

B  4  COEFFICIENTS 
89 . 68 1 5E-03  28 . 3663E-03 

ft  5  COEFFICIENTS 
-14. 1635E-03  -81 . 2760E-04 

B  5  COEFFICIENTS 
12.0177E-03  18.6757E-03 

A  6  COEFFICIENTS 
29 . 3646E-03  -98 . 1 474E-04 

B  6  COEFFICIENTS 
96 . 1 625E-04  5 1 . 9469E-04 

ft  7  COEFFICIENTS 
-56 . 6828E-04  24 . 7467E-04 

B  7  COEFFICIENTS 
-75.5487E-04  63.9953E-04 

ft  8  COEFFICIENTS 
52 . 7498E-04  -10. 3643E-04 

B  8  COEFFICIENTS 
56 . 998 1 E-04  46 . S455E-05 

ft  9  COEFFICIENTS 
-24 . 7636E-05  67 . 0423E-05 

B  9  COEFFICIENTS 
-53.9079E-04  32.3370E-04 
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010  COEFFICIENTS 


-12.1646E-04  77.5763E-05 

BIO  COEFFICIENTS 
39 . 6596E-04  1S.0515E-06 

All  COEFFICIENTS 
26 . 4455E-04  -51.1 377E-05 

81 1  COEFFICIENTS 

-42. 3032E-04  29. 1320E-04 

012  COEFFICIENTS 
-10. 2502E-03  37.9081E-04 

812  COEFFICIENTS 

38 . 4 1 64E-04  15. 5023E-04 

013  COEFFICIENTS 
53 . 7528E-04  27.8451 E-04 

813  COEFFICIENTS 

34 . 0985E-04  63 . 6847E-04 

014  COEFFICIENTS 
-36. 9928E-03  27. 0890E-03 

814  COEFFICIENTS 

29 . 8939E-03  93 . 30 1 7E-04 

015  COEFFICIENTS 
1 2 . 9997E-02  -25.201 SE-03 

815  COEFFICIENTS 


31.4443E-03 


-15. 7158E-05 


flit.  COEFFICIENTS 


10.9472E-02  -11.5703E-02 

Bit*  COEFFICIENTS 
12. 7973E-02  -14. 1906E-02 

ft 17  COEFFICIENTS 
50.5546E-02  -26.2355E-02 

B17  COEFFICIENTS 
•57 . 2088E-92  -46 . 7724E-02 

018  COEFFICIENTS 
-10. 0623E-91  45.9204E-02 

B18  COEFFICIENTS 


39.6555E-1 1 


-12. 1956E-10 


INCREMENT  OF  RUN  EQUALS  360  STEPS 
FOR  THIS  RUN; 

B=  .97  MU=  2.40  SAMMA=  6.00  PSQUARE=  1.00 
POINCARE  EXPONENTS 

OMESA 1 =- 1 5 . 49838E-0 1  0MEGA2®  86.36900E-03 

SVSTEM  OF  EIGENUECTORS 

2 1 . 2337E-02  35 . 1 497E-02 

-97.7197E-02  93.6189E-02 

A  8  COEFFICIENTS 

39 . 1 804E-0 1  -15. 8224E-0 1 

B  0  COEFFICIENTS 

.0000E+00  .0000E+00 

A  1  COEFFICIENTS 

-17. 9501E-01  68.9259E-02 

B  1  COEFFICIENTS 

-19.2140E-01  -12.9562E-01 

A  2  COEFFICIENTS 

-40. 7326E-02  49 . 9553E-02 

B  2  COEFFICIENTS 

-32.4624E-02  -31.4557E-92 

A  3  COEFFICIENTS 

-65 . 9545E-02  1 2 . 2059E-02 

B  3  COEFFICIENTS 

28. 1964E-02  47.9613E-03 


«  4  COEFFICIENTS 
22 . 0567E-02  -10.2411 E-02 

B  4  COEFFICIENTS 
14. 1571E-62  18. 4471E-03 

P  5  COEFFICIENTS 
-22.4732E-03  -16. IS52E-03 

8  5  COEFFICIENTS 
2 1 . 0586E-03  1 5 . 6389E-03 

p  6  COEFFICIENTS 
6 1 . 9723E-03  -15. 3209E-03 

B  6  COEFFICIENTS 
25 . 304 1 E-04  44 . 1 703E-04 

P  7  COEFFICIENTS 
-15.0067E-03  29.9798E-04 

B  7  COEFFICIENTS 
-15.1 755E-03  88 . 1 646E-04 

P  8  COEFFICIENTS 
10.6436E-03  -15.7823E-04 

B  8  COEFFICIENTS 
78. 9533E-04  38. 8224E-05 

ft  9  COEFFICIENTS 
-19.3759E-04  S8.2676E-05 

B  9  COEFFICIENTS 
-10. 24 04E-03  43 . 8635E -04 
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01O  COEFFICIENTS 
-22 . 2283E-04  10.1 253E-84 

BIO  COEFFICIENTS 
63. 1807E-04  -16.6157E-05 

fill  COEFFICIENTS 
56.3769E-04  -63. 4332E-05 

Bll  COEFFICIENTS 
-84 . 3872E-04  40. 1213E-04 

012  COEFFICIENTS 
-21.4691E-03  57.6457E-04 

B12  COEFFICIENTS 
22 . 4074E-04  12.01 30E-04 

013  COEFFICIENTS 
84 . 5376E-04  54.581 4E-04 

B13  COEFFICIENTS 
57 . 9542E-04  54 . 3429E-04 

014  COEFFICIENTS 
-74 . 4237E-03  34.541 0E-03 

B14  COEFFICIENTS 
47 . 4022E-03  59 . 848SE-04 

015  COEFFICIENTS 
22. 0386E-02  -40 . 7304E-O3 

B15  COEFFICIENTS 
93. 7542E-03  15. 9163E-03 
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.«rt.  «S.iA 


016  COEFFICIENTS 
13. 5244E-02  -16.6377E-02 

B16  COEFFICIENTS 
-1O.8323E-02  -1O.4925E-02 

fill?  COEFFICIENTS 
59. 8508E-02  -22.9823E-02 

B17  COEFFICIENTS 
-64 . 046 1 E-02  “43 .1931 E-02 

018  COEFFICIENTS 
-13. 0636E-0 1  52 . 745 1 E-02 

B18  COEFFICIENTS 
12. 8949E-10  -12.5211E-10 


I  m 


INCREMENT  OF  RUN  EQUALS  3b  STEPS 

B=F°?S7HIb  RmI>  2.40  SAMMA=  8.00  PSGUARE=  1.00 
POINCARE  EXPONENTS 


OMEGA 1 =-20 . 36572E-0 1 


0MEGA2=  95.5870OE-03 


SVSTEM  OF  EIGENOECTORS 


20. 5712E-02 
-97. 8613E— 02 

A  0  COEFFICIENTS 

38. 3032E-01  -83.5483E-01 

E:  0  COEFFICIENTS 


43 . 2509E-02 
90. 1630E-92 


. 0000E+00 


. 0000E+00 


19. 4091E-0S 


A  1  COEFFICIENTS 
33  41 . 5386E-01 


E:  1  COEFFICIENTS 


18. 0379E-02 


14. 8895E-01 


A  2  COEFFICIENTS 


-53. 2736E-02 


15. 6524E-01 


B  2  COEFFICIENTS 
44.7947E-01  -44.5186E-01 

A  3  COEFFICIENTS 
72 . 1 959E-0 1  -56.231 9E-0 1 

B  3  COEFFICIENTS 


-39. 0440E-01 


20. 1968E-01 
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ft  4  COEFFICIENTS 
-58 . 8626E-0 1  55 . 8466E-0 1 

B  4  COEFFICIENTS 
-2 1 . 8365E-8 1  61. 3445E-02 

ft  5  COEFFICIENTS 
17.5566E-01  -16.0705E-01 

B  5  COEFFICIENTS 
25 . 0245E-0 1  -32 . 1 868E-0 1 

ft  6  COEFFICIENTS 
46 . 7009E-0 1  -51. 0982E-0 1 

B  6  COEFFICIENTS 
-48. 7196E-01  48. 2055E-01 

ft  7  COEFFICIENTS 
-60. 263 IE -01  51 . 2633E-01 

E:  7  COEFFICIENTS 
29 . 5 1 55E-0 1  -33.31 90E-0 1 

ft  8  COEFFICIENTS 
2 1 . 1 470E-0 1  -28 . 6282E-0 1 

B  8  COEFFICIENTS 
3 1 . 588 1 E-0 1  -31. 3069E-0 1 

ft  9  COEFFICIENTS 
-32. 3638E-02  -74. 7158E-02 

B  9  COEFFICIENTS 
-45. 7714E-01  47.4797E-01 
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010  COEFFICIENTS 
-62. 3159E-01  48. 8875E-01 

B 10  COEFFICIENTS 
49 . 0673E-O i  -42 . S49SE-0 1 

Oil  COEFFICIENTS 
40 . 1 700E -0 1  -45 . 3392E-0 1 

Bll  COEFFICIENTS 
1 5 . 7928E-02  50 . 1 239E-02 

012  COEFFICIENTS 
-10.201  IE-01  67 . 6856E-02 

B12  COEFFICIENTS 
-32. 1 101E-01  42. 8889E-01 

013  COEFFICIENTS 
-26. 8316E-01  22.21 19E-01 

B13  COEFFICIENTS 
84 . 2448E-0 1  -68 . 6385E-0 1 

014  COEFFICIENTS 
36.6336E-01  -33.4362E-01 

B14  COEFFICIENTS 
-11. 6566E-0 1  19. 4679E-0 1 

015  COEFFICIENTS 
-42 . 9 1 79E-0 1  43.01 88E-0 1 

B15  COEFFICIENTS 


INCREMENT  OF  RUN  EQUALS  360  STEPS 
FOR  THIS  RUN; 

B=  .97  MU=  2.40  SAMMA=10.00  PSQUARE=  1.00 

POINCARE  EXPONENTS 

OMEGA 1 =-25. 12578E-01  0MEGA2=  81.90400E-03 

SVSTEM  OF  EIGENUECTORS 

-18.31 84E-02  50 . 7958E-02 

98. 3079E-02  86. 1382E-02 

A  0  COEFFICIENTS 

78. 1317E-01  -18. 7741E-01 

G  0  COEFFICIENTS 

.0000E+00  .0000E+00 

A  1  COEFFICIENTS 

35.9970E-01  54.7646E-02 

B  1  COEFFICIENTS 

34. 0609E-01  -14.3619E-01 

A  2  COEFFICIENTS 

27 . 2329E-02  67 . 4262E-02 

B  2  COEFFICIENTS 

38 . 9324E-02  -20 . 0546E-02 

A  3  COEFFICIENTS 

17. 2081E-01  15.4650E-02 

B  3  COEFFICIENTS 

•10.6336E-01  15.0248E-82 


•W 


4  ■j  M  M  ■  ■  ■  1 


ft  4  COEFFICIENTS 
83. 4205E-02  -14. 0620E-02 

B  4  COEFFICIENTS 
30 . 4409E-02  17.1 890E-83 

ft  5  COEFFICIENTS 
1 0 . 4539E-02  -28 . 8765E-03 

B  5  COEFFICIENTS 
78.2825E-03  89. 7957E-04 

ft  6  COEFFICIENTS 
-22. 6880E-02  -25. 1819E-03 

B  6  COEFFICIENTS 
70. 281 0E-03  56 . 6429E-04 

ft  7  COEFFICIENTS 
84 . 7070E-03  27. 9216E-04 

B  7  COEFFICIENTS 
39 . 0608E-03  1 6 . 2973E-03 

Pi  8  COEFFICIENTS 
-38 . 30 1 6E-03  -33 . 7700E-04 

B  8  COEFFICIENTS 
-14.51 42E-03  35 . 6963E-04 

PI  9  COEFFICIENTS 
15. 6598E-03  86.7653E-05 

B  9  COEFFICIENTS 
25. 5472E-03  89.6789E-04 
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Pile  COEFFICIENTS 


52 . 2396E-04  92.451 8E-05 

810  COEFFICIENTS 

-17.6781E-03  20.3526E-84 

Pill  COEFFICIENTS 
-25 . 9063E-03  -10. 420 1 E-04 

811  COEFFICIENTS 

23 . 2306E-03  8 1 . 9764E-04 

PI  12  COEFFICIENTS 
77 . 0499E-03  85 . 9993E-94 

812  COEFFICIENTS 

1 6. 7042E-03  24.1281 E-94 

PI  13  COEFFICIENTS 
-36 . 6687E-03  9 1 . 0647E-04 

813  COEFFICIENTS 

-22 . 0450E-03  4 1 . 8428E-04 

PI14  COEFFICIENTS 
23.0792E-02  46.8645E-03 

814  COEFFICIENTS 

-10.3866E-02  61.0813E-04 

PI  15  COEFFICIENTS 
-57 . 536 1 E-02  -52 . 1 643E-03 

815  COEFFICIENTS 

-35 . 3607E-02  50 . 5087E-03 


N16  COEFFICIENTS 


88.9147E-03  -22.5051E-02 

616  COEFFICIENTS 
12.9291E-02  -66.6425E-03 

PI17  COEFFICIENTS 
■12.0094E-01  -18.3170E-02 

B17  COEFFICIENTS 
11.3537E-01  -47 . 8633E-02 

fil8  COEFFICIENTS 
26 . 0573E-6 1  62 . 5389E-02 

618  COEFFICIENTS 
-47.2771E-10  -13.1795E-10 


INCREMENT  OF  RUN  EQUALS  360  STEPS 

FOR  THIS  RUN;  n„„.DC  ,  ^ 

B=  .97  MU=  2.40  6AMMA=  1 2 . Ob'  PSQUHRE*  1.0* 

POINCARE  EXPONENTS 

OMEGA 1 =-29. 86568E-01  0MEGA2=  56.85900E-03 

SVSTEM  OF  EIGENUECTORS 

-15.1443E-02  5? . 8038E-02 

98 . 8466E-02  8 1 . 60 1 OE -02 

A  O  COEFFICIENTS 

-11.0410E+O0  -29.2664E-01 

B  O  COEFFICIENTS 

. 00O0E+00  ■ 0000E+0# 

A  1  COEFFICIENTS 

53.61 73E-0 1  46 . 0498E-02 

B  1  COEFFICIENTS 

44 . 5590E-0 1  -30 . 6984E-02 

A  2  COEFFICIENTS 

-20.3946E-02  20.7150E-01 

B  2  COEFFICIENTS 

55 . 2562E-02  79 . 6890E-02 

A  3  COEFFICIENTS 

26.81 1 IE-01  18.6141E-01 

B  3  COEFFICIENTS 


-17.2685E-01 


76. 1854E-03 


ft  4  COEFFICIENTS 
-14.8115E-01  10. 3161E-01 

B  4  COEFFICIENTS 
-32 . 39 1 3E-62  -83 . 7459E-02 

ft  5  COEFFICIENTS 
24 . 2456E-02  54 . 4349E-02 

B  5  COEFFICIENTS 
-15. 3655E-02  -95 . 0396E-02 

ft  6  COEFFICIENTS 
-39. 6881E-02  20. 8239E-02 

B  6  COEFFICIENTS 
13. 1990E-02  -76. 2580E-02 

ft  7  COEFFICIENTS 
1 7 . 4706E-02  1 2 . 3798E-02 

B  7  COEFFICIENTS 
46 . 6590E-03  -55 . 9 1 42E-02 

ft  8  COEFFICIENTS 
-71.121 9E-03  90 . 6268E-03 

B  8  COEFFICIENTS 
-95.5110E-04  “46. 1019E-02 

ft  9  COEFFICIENTS 
33.6707E-03  37.6234E-03 

B  9  COEFFICIENTS 
35.4174E-03  -39.6964E-02 
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me 

COEFFICIENTS 

61.3447E-04 

78. 7251E-03 

Bie 

COEFFICIENTS 

22.3023E-03 

-37.9507E-02 

mi 

COEFFICIENTS 

50. 9039E-03 

60. 1988E-03 

Bll 

COEFFICIENTS 

33.4580E-03 

-37. 3592E-02 

m2 

COEFFICIENTS 

13. 2496E-02 

25.0669E-83 

B12 

COEFFICIENTS 

49.5796E-03 

-41.2900E-02 

013 

COEFFICIENTS 

82. 8849E-03 

-94.0776E-63 

B13 

COEFFICIENTS 

42 . 3968E-03 

-44.4325E-02 

014 

COEFFICIENTS 

49. 7726E-02 

-26.0973E-92 

B14 

COEFFICIENTS 

11. 2193E-02 

-38.0065E-02 

015 

COEFFICIENTS 

89.6710E-02 

-54.2180E-02 

B15 

COEFFICIENTS 

■57.3149E-02 

-49. 1 236E-03 

G16  COEFFICIENTS 


71 . 1015E-03  -61.4583E-02 

B16  COEFFICIENTS 
18.3251 E-02  2 1 . 6284E-02 

B17  COEFFICIENTS 
■17.8941E-01  -79.61S3E-03 

B17  COEFFICIENTS 
14.8565E-01  -12.6878E-02 

fll8  COEFFICIENTS 
36.8272E-01  10.4914E-01 

B18  COEFFICIENTS 
-79 . 5457E- 10  -11. 3262E- 1 0 


INCREMENT  OF  RUN  EQUALS  368  STEPS 
FOR  THIS  RUN; 

B=  .97  MU=  2.48  GAMMA=14.00  PSQUARE=  1.00 

POINCARE  EXPONENTS 

OMEGAl=-30. 31760E-01  0MEGA2=  25.30000E-03 

SVSTEM  OF  EIGENUECTORS 

-11. 3653E-02  64. 3938E-02 

99.3521E-02  76.5077E-02 

A  0  COEFFICIENTS 

-85 . 9898E+00  -19.141 0E+08 

B  0  COEFFICIENTS 

. 0000E+00  - 0000E+00 

A  1  COEFFICIENTS 

26.81 25E+80  34 . 9980E-0 1 

B  1  COEFFICIENTS 

48 . 485 1 E+00  27 . 8354E+00 

A  2  COEFFICIENTS 

1 2. 5882E+00  37. 0394E+00 

B  2  COEFFICIENTS 

68.6565E-01  23.7821E+00 

A  3  COEFFICIENTS 

24. 7123E+00  47.2115E+80 

B  3  COEFFICIENTS 

-16.3433E+00  -34. 4419E-01 


ft  4  COEFFICIENTS 
-11.0840E+00  34.4015E+00 

B  4  COEFFICIENTS 
-97.9130E-01  -24.2313E+00 

ft  5  COEFFICIENTS 
-70 . 0527E-02  1 7 . 4554E+00 

B  5  COEFFICIENTS 
-22. 7835E-0 1  -28. 1987E+00 

ft  6  COEFFICIENTS 
-41 . 3841E-01  70. 2036E-01 

B  6  COEFFICIENTS 
98. 6141E-02  -23.2191E+00 

ft  7  COEFFICIENTS 
12. 0444E-01  32.4451E-01 

B  7  COEFFICIENTS 
50.0426E-02  -17.4411E+08 

ft  8  COEFFICIENTS 
-48 . 342 1 E-02  24 . 6653E-0 1 

B  8  COEFFICIENTS 
-40. 7154E-02  -13.8133E+00 

ft  9  COEFFICIENTS 
29.3130E-02  24.1659E-01 

B  9  COEFFICIENTS 
-71.3964E-03  -11.9622E+00 


010  COEFFICIENTS 
31.171 4E-03  23 . 2357E -9 1 

B10  COEFFICIENTS 
-5 1 . 3457E-02  -11. 2468E+00 

fill  COEFFICIENTS 
-34 . 7540E-02  1 8 . 9489E-0 1 

611  COEFFICIENTS 

1 1 . 593 1  E-02  -11. 4030E+00 

012  COEFFICIENTS 
13. 7792E-0 1  44. 3689E-02 

612  COEFFICIENTS 

89. 1242E-Q3  -12. 4071 E +00 

013  COEFFICIENTS 
2 1 . 3890E-02  -32 . 1 465E-0 1 

613  COEFFICIENTS 

-82. 2281 E-02  -13.2095E+00 

014  COEFFICIENTS 
37 . 1 995E-0 1  -90.111 2E-0 1 

614  COEFFICIENTS 

-33. 9488E-01  -11.0791E+00 

015  COEFFICIENTS 
-82.6533E-01  -13.3986E+99 

615  COEFFICIENTS 

-55.0265E-01  -33. 7420E-01 


PI  16  COEFFICIENTS 


■41 . 7859E-01  -IS. 0905E+00 

B16  COEFFICIENTS 
22 . 3279E-0 1  64 . 5452E-0 1 

Pi  17  COEFFICIENTS 
■89. 6017E-01  10. 3887E-01 

B17  COEFFICIENTS 
16.1 402E+00  85 . 4520E-0 1 

ft 18  COEFFICIENTS 
28 . 6747E+00  85 . 6943E-0 1 

B18  COEFFICIENTS 
■67 . 6235E-09  3S.6049E-10 


INCREMENT  OF  RUN  EQUALS  360  STEPS 
FOR  THIS  RUN; 

.97  MU=  2.40  GAMMA=16.00  PSQUARE=  1.00 

POINCARE  EXPONENTS 

OMEGA 1 =-33. 36009E-01  OMEGA2=-10. 56400E-03 

SVSTEM  OF  EIGENUECTORS 

-70. 7829E-03  70.5935E-02 

99. 7492E-02  70. 8276E-02 

A  0  COEFFICIENTS 

-23. 8665E+01  -13.7204E+O1 

B  0  COEFFICIENTS 

. 0000E+00  . 0000E+00 

A  1  COEFFICIENTS 

64. 4310E+09  73.7865 E+00 

B  1  COEFFICIENTS 

14.1 890E+0 1  31. 0854E+0 1 

A  2  COEFFICIENTS 

42.8193E+00  43.5389E+01 

B  2  COEFFICIENTS 

25. 1 191E+00  26.0912E+01 

A  3  COEFFICIENTS 

75. 7699E+00  56. 9025E+01 

B  3  COEFFICIENTS 

-49 . 0324E+00  -45 . 6522E+00 


ft  4  COEFFICIENTS 


-32 . 8602E+00  43 . 5564E+0 1 

B  4  COEFFICIENTS 
-33. 1830E+00  -29. 5376E+01 

ft  5  COEFFICIENTS 
-44.8871E-01  22. 9534E+01 

B  5  COEFFICIENTS 
-76. 0662E-01  -35. 7071E+81 

ft  6  COEFFICIENTS 
-13. 9459E+60  92 . 3757E+00 

B  6  COEFFICIENTS 
4 1 . 4565E-0 1  -30 • 1 447E+0 1 

ft  7  COEFFICIENTS 
38.2142E-01  33.4888E+00 

B  7  COEFFICIENTS 
22. 6164E-01  -22.6983E+01 

ft  8  COEFFICIENTS 
-13.4716E-01  27.8865E+00 

B  8  COEFFICIENTS 
-75.0055E-02  -17.6S13E+01 

ft  9  COEFFICIENTS 
1 0 . 5338E-0 1  29.141 5E+00 

B  9  COEFFICIENTS 
-24.2361E-02  -15. 1420E+01 


COEFFICIENTS 


84.3313E-04 

810 

■13. 3666E-0 1 

fill 

■10.  75&1E-01 

81 1 

60. 809 IE-02 

ft  1 2 

45.9554E-01 

812 

68. 7915E-02 

R13 

14. 1300E-01 

813 

■26. 9712E-01 

014 

11.0020E+00 

814 

■11. 4450E+00 

ftl5 

■25.3529E+0O 

815 

-16.4896E+00 


29.9560E+0O 
COEFFICIENTS 
-14. 2526E+01 
COEFFICIENTS 
25. 0O97E+0O 
COEFFICIENTS 
-14. 6181E+01 
COEFFICIENTS 
48.6618E-01 
COEFFICIENTS 
-15.9560E+01 
COEFFICIENTS 
-43.0647E+00 
COEFFICIENTS 
-16. 7415E+01 
COEFFICIENTS 
-11.3643E+01 
COEFFICIENTS 
-13.6665E+01 
COEFFICIENTS 
-15.9643E+01 
COEFFICIENTS 
-43.5968E+00 


016  COEFFICIENTS 


-14.2387E+00  -11.6184E+01 

B16  COEFFICIENTS 
82 . 0835E-0 1  68.1 757E+00 

017  COEFFICIENTS 
-2 1 . 5640E+00  36.931 7E-0 1 

E17  COEFFICIENTS 
47 . 2386E+08  94 • 2555E+00 

018  COEFFICIENTS 
79. 5737E+00  73. 8045E+00 

818  COEFFICIENTS 
-20. 4307E-08  58.6715E-09 
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100  REM  4TH  ORDER  RUNGA-KUTTA  ROUTINE  FOR  UP  TO  AN  8TH  ORDER  SVSTEM 
11©  REM 

120  CALL  CLEAR 
130  REM 

140  OPEN  #1: "PIG" 

150  REM 

160  DIM  X<8, 8) ,  XW0RKC8, 8) ,  XNEU<8, 8)  ,  AMAT<S“,  8> 

170  DIM  Fl<3.8) ,F2<8.8>,F3<8,8> ,F4C8.8>,F<8,8> 

180  DIM  K1 (8, 8>,K2<8,8>,K3< 8, 8>,K4<8. 8) 

190  REM 

200  REM  INITIALIZE  ALL  UARIABLES 

210  REM 

220  AN6LE=0.0 

230  FOR  1=1  TO  8 

240  FOR  J=1  TO  8 

250  X<I,J>=0.0 

260  IF  I*J  THEN  X<I,J>*1.8 

270  NEXT  J 

280  NEXT  I 

290  REM 

300  CALL  CLEAR 
310  REM 

320  REM  ENTER  THE  INPUT  PARAMETERS 
330  REM 

340  INPUT  "ENTER  «  OF  INCREMENTS": I NCR 
350  DELTAN6LE=2*P I / I NCR 
360  REM 

370  INPUT  "ENTER  ORDER  OF  SVSTEM": ORDER 

372  INPUT  "ENTER  IN  B:":BE 

374  INPUT  "ENTER  IN  MU:": MU 

376  REM  INPUT  "ENTER  IN  GAMMA: GAMMA 

378  INPUT  “ENTER  IN  PSQUARE: “ : PSQUARE 

379  CALL  CLEAR 

380  REM 

390  REM  MASTER  LOOP 
400  REM 

405  FOR  6AMMA=2  TO  16  STEP  2 

406  DISPLAV  AT<5,4>:"GAMMA=\ GAMMA 
410  FOR  N=1  TO  I NCR  STEP  1 

420  DISPLAV  AT<2,4): "N=" ,N 
430  REM 

440  REM  FIRST  RUNGA-KUTTA  CONSTANT 
450  REN 

460  CALL  FCNN < ANGLE ,X(,),FK,  >„BE,MU, GAMMA , PSQUARE ) 

470  REM 

480  FOR  1=1  TO  ORDER 
490  FOR  J=1  TO  ORDER 
500  K1CI, J>=DELTANGLE*F1<I, J) 

530  REM 


566  XUIORKC  I . J>=XC  I ,  J)+e.  5*K  1  <  I .  J> 

570  NEXT  J 
580  NEXT  I 
596  REM 

600  REM  SECOND  RUNGA-KUTTA  CONSTANT 
610  REM 

620  CALL  FCNN C ANGLE +0.5*DEL TANGLE, XWORKC, >,F2C.  >, BE, MU, GAMMA, PSQUARE) 
630  REM 

640  FOR  1=1  TO  ORDER 
650  FOR  J=1  TO  ORDER 
660  K2  < I , J ) =OELT AN6LE*F2  <  I ,  J  > 

69©  REM 

720  XWORK < I , J > =X <1 , J >  +0 . 5*K2 <  I ,  J  > 

730  NEXT  J 
740  NEXT  I 
750  REM 

760  REM  THIRD  RUNGA-KUTTA  CONSTANT 
770  REM 

780  CALL  FCNN C AN6LE+0 . 5*DELT ANGLE , XUORK < , > , F3C , > , BE, MU, GAMMA , PSQUARE ) 
790  REM 

800  FOR  1=1  TO  ORDER 
810  FOR  J=1  TO  ORDER 
820  K3CI, J>=DELTANGLE*F3CI. J) 

850  REM 

880  XUIORKC  I,  J)=XCI, J)+K3CI, J> 

890  NEXT  J 
900  NEXT  I 
910  REM 

920  REM  FOURTH  RUNGA-KUTTA  CONSTANT 
930  REM 

940  CALL  FCNN  C  ANGLE+DELTANGLE , XWORK  C , ) , F4  < , ) , BE , MU ,  GAMMA , PSQUARE  > 

950  REM 

960  FOR  1=1  TO  ORDER 
970  FOR  J=1  TO  ORDER 
980  K4  < I , J ) =DELTANGLE*F4  < I , J ) 

990  NEXT  J 
1000  NEXT  I 
1010  REM 

1020  REM  EUALUATION  OF  STATE  VARIABLES 
1038  REM 

1040  FOR  1=1  TO  ORDER 
1050  FOR  J=1  TO  ORDER 

1060  XNEWC I , J>=X< I , J)+< 1/6)*<K1 < I , J)+2*K2< I , J>+2*K3< I , J>+K4< I , J) ) 

1070  NEXT  J 

1080  NEXT  I 

1090  REM 

1100  REM  UPDATE  STATE  VARIABLES 
1110  REM 

1120  FOR  1=1  TO  ORDER 
1130  FOR  J=1  TO  ORDER 
1140  X  < I , J  > =XNEU  < I , J  > 

1150  NEXT  J 

1160  NEXT  I 
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1 360  SUE  FCNN  <  T ,  X WORK  C ,  > ,  F  C ,  ) ,  BE ,  MU ,  6AMMA ,  PSQUARE  > 

1370  REM 

1 3S0  CALL  MATR I XA  <  T , AMAT  < . > , BE , MU . GAMMA , PSQUARE > 

1390  REM 

1400  CALL  MATRIXMULTCAMATC,  >,XUORKC,  >,2,2,2,FC, >) 

1410  REM 
1420  REM 
•1430  SU8END 

2000  SUB  MATRIXACPSI , MAT AC , >,BE, MU, GAMMA, PSQUARE > 

2001  REM 

2002  REM  ALGORITHM  TO  CALCULATE  STATE  TRANSITION  MATRIX 

2003  REM  FOR  A  PERIODIC  FUNCTION 

2004  REM  IN  PARTICULAR-HELICOPTER  BLADE,  NO  FEEDBACK 

2005  REM  UPDATED  16  AUG  1984 

2006  REM 

2007  DEF  ARCSINCZZ>=ATNCZZ/SGRC1-ZZ*ZZ>> 

2008  PSIT=PSI 

2010  IF  PSIT>=0.0  AND  PS I TCP I  THEN  2020  ELSE  2012 

2012  IF  PSIT>=PI  AND  PS I TC < P I +ARCS I N ( BE^MU > ) THEN  2030  ELSE  2014 

2014  IF  PS  I  T>=  <  2*P  I  -ARCS  I NC  BE/MU  )  >  AND  PSIT«2*PI>-.  008001  THEN  283©  ELSE  2016 

2016  IF  PS I T>C  2*P I ) 00900 1  THEN  PSIT=0.0  ::  GO  TO  2020  ELSE  60  TO  2018 

2018  GO  TO  2040 

2020  REM  PRINT  « Is" 2020" 

2022  K=  <1 /Z ) *6E~3*MU*C0S C  PS I T ) + C 1 ^4  > *BE A2*MUA2*S I N < 2*PS I T > 

2023  REM 

2024  C=  Cl /4 > *6EA4+  C 1 /Z ) *6EA3*MU*S I N C PS I T ) 

2025  REM 

2026  GO  TO  2050 

2030  REM  PRINT  #1:"2030" 

2031  REM 

2032  K=Cl/3>*B£'C5*MU*C0SCPSm-Kl/^4>*BE''2*MU''2*SINC2*PSIT>+MLr4*C-Cl''12>*SIN<2*P 
SIT)+C 1/24)*SINC4*PSIT) ) 

2033  REM 

2034  C=C1/4>*BEA4+C1/3>*BE''3*MU*SINCPSIT)+MUA4*CC1/'16)-C1^12>*C0SC2*PSIT)+C1/4B> 
*C0SC4*PSIT)> 

2035  REM 

2036  GO  TO  2050 

2040  REM  PRINT  # Is" 2040“ 

2041  REM 

2042  K=- C 1  /3 > *BE/'3*MU*CQS CPSIT)-Cl/4) *BE/'2*MU/'2*S I N C 2*PS I T > 

2043  REM 

2044  C=-  C 1  /4  )  *BE  A4-  C 1  ^3  >  *BE'‘3*MU*S  I N  <  PS  I T  ) 

2045  REM 

2050  REM  CONSTRUCT  MATRIX  A 

2051  REM 

2052  MATAC1, 1>=0 

2053  MATAC 1 , 2>=1 

2054  MATAC2, 1)»-C6AMMA/-2>*C2*PSQUARE/GAMMA+K> 

2055  MATAC2,2)=-CGAMMA*C/2> 

2056  REM 

2057  REM  RETURN  TO  MAIN  PROGRAM 
2060  SUBEND 
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1170  REM 

1180  ANGLE= ANGLE +OEL  T ANGLE 
1190  REM 
1200  NEXT  N 
1210  REM 

1220  REM  PRINT  OUTPUT 

1230  REM 

1231  PRINT  #1:“ - 

11 

1232  PRINT  #1 

1233  PRINT  #1, USING  "  BE=##.  MU=##.  :  BE, MU 

1234  PRINT  #1 

1235  PRINT  #1, USING  “  ,GAMMA=##.  PSGlUARE=##.  : GAMMA,  PSQUARE 

1236  PRINT  #1 

1237  PRINT  #1 

1240  FOR  1=1  TO  ORDER 
1250  FOR  J=1  TO  ORDER 

1260  PRINT  #1, USING  “  X <#,#>=###. : I, J,X(I, J> 

1270  NEXT  J 
1280  NEXT  I 
1290  PRINT  #1 
1300  PRINT  #1 
1310  PRINT  #1 
1320  REM 

1325  NEXT  GAMMA 

1326  REM 
1330  CLOSE  #1 
1340  REM 
1350  END 


•.  - v\% 
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7000  SUE'  MATR  1 XMUL T(£K,  >,B<,  >,  ROUIM . COLA .  COLB ,  OUT < .  > > 

7801  REM 

7002  REM  MATRIX  MULTIPLICATION  ROUTINE 

7003  REM  INPUTS  MATRIX  “A".  MATRIX  ■'B" . DIMENSION  OF  AAND  COLUMN  DIMENSION  OF 
B 

7004  FOR  1=1  TO  ROWA 

7005  FOR  J=1  TO  COLA 

7006  OUT < I, J >=0.0 

7007  NEXT  J 

7008  NEXT  I 

7009  REM 

7010  FOR  1=1  TO  ROWA 
7020  FOR  K=1  TO  COLB 
7030  FOR  J=1  TO  COLA 
7040  SUM=A< I , J)*B< J, K> 

7050  OUT < I , K)=GUT < I , K>+SUM 
7060  NEXT  J 

7070  NEXT  K 
7080  NEXT  I 
7090  REM 

7999  SUBEND  ::  REM  RETURN  TO  MAIN  PROGRAM 
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100  REM  PROGRAM  EIGENP-  EIGENUALUE/EIGENUECTOR  ROUTINE 
110  REM 

120  OPEN  #1: "PIO" 

125  CALL  CLEAR 

130  DIM  A<4, 4) , UECRC4, 4>, UECI <4, 4), EUR<4), EUI <4> , INDIC<4) 

140  DIM  I WORK ( 10) , LOCAL< 10) , PRFACT < 10) , SUBDI A( 10) 

150  DIM  UORKl(10),WORK2<10),ldORK<10),H<‘4,4) 

170  REM 

180  PRINT  “ENTER  IN  A  MATRIX" 

190  CALL  INPUTT(N,NM,A<, ), IERROR) 

200  REM 
210  T=64 
220  REM 

230  FOR  1=1  TO  N 
235  FOR  J=1  TO  N 

240  PRINT  #1,  USING  "  A<##,  ##)=##.  " :  I ,  J,  A< I ,  J) 

245  NEXT  J 
250  NEXT  I 

255  REM 

256  PRINT  #1 

257  PRINT  #1 

258  PRINT  #1 
260  REM 

360  REM 

370  IF  NOl  THEN  450 
380  EUR(l)=Aa,l) 

390  EUI < 1 )=0. 0 

400  UECR<1,1)=1.0 

410  UECI<1,1>=0.0 

420  IND1C<1)=2 

430  GO  TO  2080 

440  REM  STATEMENT  LABEL  1 

450  CALL  SCALEE<N,NM,A<, ),UECI<, ), PRFACT O , ENORM) 

460  REM  THE  COMPUTATION  OF  THE  EIGENUALUES  OF  A  NORMALIZED  MATRIX 
470  REM 

480  EX=EXP<-T*LOG<2.0)> 

490  CALL  HESQR<N, NM, A< , ) , UECI < , ) , EUR< ) , EUI O , SUBDIA< ) , INDIC< ) , EPS, EX) 
500  J=N 
510  1=1 

520  LOCAL< 1 )=1 

530  IF  J=1  THEN  650 

540  REM  STATEMENT  LABEL  2 

550  IF  ABS(SUBDIA< J-l >>>EPS  THEN  580 

560  1=1+1 

570  LOCAL < I )=0. 

580  REM  STATEMENT  LABEL  3 
590  J=J-1 

600  LOCAL  < I ) =LOC AL  < I ) + 1 
610  IF  JOl  THEN  540 


620  REM 

630  REM  EIGENVECTOR  PROBLEM  / 

640  REM 

650  REM  STATEMENT  LABEL  4 
660  K=1 
670  KGN*0. 

680  L*LOCAL< 1 > 

690  M=N 

708  FOR  1=1  TO  N 

710  IUEC=N-I+1 

720  IF  IOL  THEN  760  ' 

73©  K=K+1  i 

740  M=N-L  j 

750  L=L+LOCAL<K)  , 

760  REM  STATEMENT  LABEL  5  I 

770  IF  INDIC<IUEC)=8  THEN  1070 
780  IF  EUIUUEOO0.0  THEN  1000 
790  REM 

800  REM  TRANSFER  OF  AN  UPPER-HESSENBER6  MATRIX  OF  THE  ORDER 
810  REM  M  FROM  THE  ARRAVS  VECI  ANO  SUBDIA  INTO  ARRAV  A 
820  REM 

830  FOR  Kl=l  TO  M 

840  FOR  H=K1  TO  M 

850  A<K1,L1)=UECI<K1,L1) 

860  NEXT  LI 

870  IF  Kl=l  THEN  890 

880  A<K1,K1-1)=SUBDIA<K1-1) 

890  REM  STATEMENT  LABEL  7 
900  NEXT  K1 
910  REM 

920  REM  THE  COMPUTATION  OF  THE  REAL  EIGENVECTOR  IUEC  FOR  THE  UPPER  HESSENBER6  MA 
TRIX 

930  REM  CORRESPONDING  TO  THE  REAL  EIGENVALUE  EURCIVEC) 

940  REM 

950  CALL  REAL VE<N,NM,M,  IUEC, AC,  ),VECRC,  ), EURO, EVIO,  IWORKO,WORKO,  INDICO, EPS, 

EX) 

960  60  TO  1070 
970  REM 

980  REM  THE  COMPUTATION  OF  THE  COMPLEX  EIGENVECTOR 
990  REM 

1000  REM  STATEMENT  LABEL  8 
1010  IF  KONO0  THEN  1050 
1020  K0N=1 

1030  CALL  COMPVE<N,  NM,  M,  IUEC,  A< ,  ) ,  UECRC ,  ) ,  VECI  < ,  > ,  EVR<  ) ,  EVI  <  ) .  INDIC-C  ) ,  I  WORK  C  ) ,  SU 
BDIA<  ) ,  W0RK1  <  ) ,  W0RK2<  ) ,  WORKO ,  EPS,  EX) 

1040  GO  TO  1070 

1050  REM  STATEMENT  LABEL  9 

1060  KON=0 

1070  REM  STATEMENT  LABEL  10 
1875  NEXT  I 
1088  REM 
1090  REM 

1100  REM  THE  RECONSTRUCTION  OF  THE  MATRIX 
1110  REM 
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1 

I 


1120  FOR  1=1  TO  N 
1130  FOR  J=1  TO  N 
1140  A<I,J>=0.0 
1 150  A<  J, I >=0.0 
1160  NEXT  J 
1180  A<I,  I)=l. 

1190  NEXT  I 

1200  IF  N<=2  THEN  1370 

1210  M=N-2 

1220  FOR  K=1  TO  M 

1230  L=K+1 

1240  FOR  J=2  TO  N 

1250  D1=0.0 

1260  FOR  I=L  TO  N 

1270  D2=UECI<I,K> 

1280  D1=D1+D2*A<J,I> 

1290  NEXT  I 

1300  FOR  I=L  TO  N 

1310  A<J,I>=A<J,I>-UECI<I,K>*D1 

1320  NEXT  I 

1325  NEXT  J 

1330  NEXT  K 

1340  REM 

1350  REM  COMPUTE  THE  EIGENUECTORS  OF  THE  ORIGINAL  NON-SCALED  MATRIX 
1360  REM 

1370  REM  STATEMENT  LABEL  15 
1380  KON=l 
1390  FOR  1=1  TO  N 
1400  L=8. 

1410  IF  EUI<I)=0.0  THEN  1460 
1420  L=1 

1430  IF  KON=0. 0  THEN  1460 

1440  KQN=0 

1450  GO  TO  2010 

1460  REM  STATEMENT  LABEL  16 

1470  FOR  J=1  TO  N 

1480  01=8.0 

1490  D2=0.0 

1500  FOR  K=1  TO  N 

1510  D3=A<J,K> 

1520  D1=D1+D3*UECR<K, I> 

1530  IF  L=0  THEN  1550 
1 540  D2=D2+D3*UECR <  K . I - 1 > 

1550  NEXT  K  ::  REM  STATEMENT  LABEL  17 
1560  WORK  <  J  )  =0 1  <""PRFACT  <  J  > 

1570  IF  L=0.0  THEN  1590 
1580  SU8D I A  <  J  ) =D2/PRFACT  <  J ) 

1590  REM  STATEMENT  LABEL  18 
1600  NEXT  J 
1610  REM 

1620  REM  THE  NORMALIZATION  OF  THE  EIGENUECTORS  AND  THE  COMPUTATION  OF  THE 
ALUES 

1630  REM  OF  THE  ORIGINAL  NON-NORMAL I ZED  MATRIX. 


EIGEN' 


1640  IF  L=l  THEN  1760 
1650  D1=0.0 
1660  FOR  M=1  TO  N 
1670  D1=01+W0RK(M>''2 
1680  NEXT  M  ’ 

1696  D1=SQR<D1> 

1790  FOR  M=1  TO  N 
1710  UECKM,  I  >=0.0 
1720  OECR<M,  n=W0RK<M>/01 
1730  NEXT  M 

1740  EUR  <1 > =EUR < I > *ENORM 
1750  GO  TO  2810 
1760  REM  SRSTATEMENT  LABEL  21 
1770  KON=l 

1780  EUR< I )=EUR< I >*ENORM 
1790  EUR( 1-1 >=EUR< I > 

1800  EUI < I >=EUI ( I )*ENORM 
1810  EUI < 1-1 )=-EUI < I ) 

1820  R=0.0 

1830  FOR  J=1  TO  N 

1840  Rl=WORK(  J)/'2+SUBDIA<  J)A2 

1850  IF  R>=R1  THEN  1880 

I860  R=R1 

1870  L=J 

1880  REM  STATEMENT  LABEL  22 
1890  NEXT  J 

1900  D3=W0RKCL)  ' 

1918  R1=SUBDIA<L) 

1920  FOR  J*1  TO  N 
1930  Dl=WORK< J) 

1940  D2=SUBDIA< J) 

1950  UECRCJ, I)=<D1*D3+D2*R1VR 
I960  UECKJ,  I)=<D2*D3-D1*R1)/R 
1970  UECR<J, I-1>=UECR<J, I) 

1980  REM  STATEMENT  23 
1990  UECKJ,  I-1)=-UECKJ,  I) 

2000  NEXT  J 

2010  REM  STATEMENT  24 

2020  NEXT  I 

2070  REM 

2080  REM  STATEMENT  25 

2081  FOR  1=1  TO  N 

2082  PRINT  #1, USING  "EIGENUALUE  ##  EQUALS  +  ##. I,EUR<I>,E 

UI(I> 

2083  PRINT  #1, USING  "EI6ENUECT0R  EQUALS" 

2084  FOR  J=1  TO  N 

2085  PRINT  #1, USING  "  ### .  +  ###. #####^^JM:UECR<J,  I), UECKJ,  I> 

2086  NEXT  J 

2087  PRINT  #1 

2088  NEXT  I 

2089  CLOSE  #1 

2090  END 
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3800  SUB  SCALEE<N,NM,A<,  >,H<,  ),PRFACTO,ENORM> 
3010  REM 
3020  REM 

3930  FOR  11=1  TO  N 
3040  FOR  JJ=1  TO  N 
3050  H<II. JJ)=A<II,JJ> 

3060  NEXT  JJ 
3070  PRFACT  < 1 1 >=1 . 0 
3080  NEXT  II 
3090  BOUND 1=0. 75 
3100  B0UND2= 1 . 33 
3110  ITER=0 

3120  REM  STATEMENT  LABEL  3 
3130  NCOUNT=0. 

3140  FOR  11=1  TO  N 

3150  COLUMN=0. 0 

3160  ROW=0. 0 

3170  FOR  JJ=1  TO  N 

3180  IF  II=JJ  THEN  3210 

3 1 98  COLUMN=COLUMN+ABS (flCJJ,  ID) 

3200  ROW=ROW+ ABS <  A  <  1 1 ,  J  J  )  ) 

3210  REM  STATEMENT  4 
3220  NEXT  JJ 

3230  IF  COLUMN=0.0  THEN  3280 

3240  IF  ROW=0.0  THEN  3280 

3250  Q=COLUMN/'ROW 

3260  IF  GKBOUNDl  THEN  3310 

3270  IF  Q>BOUND2  THEN  3310 

3280  REM  STATEMENT  LABEL  5 

3298  NCOUNT=NCOUNT  + 1 

3300  60  TO  3400 

3310  REM  STATEMENT  LABEL  6 

3328  FACTOR=SQR(Q> 

3330  FOR  JJ=1  TO  N 
3340  IF  II=JJ  THEN  3370 
3350  Ad  I,  JJ)=A<II,JJ)*FACTOR 
3360  A( JJ, 1 1 )=A< JJ, 1 1 ) /FACTOR 
3370  REM  STATEMENT  LABEL  7 
3380  NEXT  JJ 

3390  PRFACT < 1 1 )=PRFACT< 1 1 >*F ACTOR 

3400  REM  STATEMENT  LABEL  8 

3410  NEXT  II 

3420  ITER=ITER+1 

3430  IF  ITER>30  THEN  3610 

3440  IF  NCOUNT<N  THEN  3120 

3450  FNORM=0.0 

3460  FOR  11=1  TO  N 

3470  FOR  JJ=1  TO  N 

3480  Q=A<II, JJ) 

3490  FNORM=FNORM+Q*Q 
3500  REM  STATEMENT  LABEL  9 
3510  NEXT  JJ 


3525  FNORM=SQR ( FNORM ) 

3530  FOR  I  I=*l  TO  N 

3540  FOR  JJ=1  TO  N 

3550  A<  1 1 ,  JJ>=A<  1 1 .  JJ  VFNGRM 

3560  REM  STATEMNET  LABEL  10 

3570  NEXT  JJ 

3580  NEXT  II 

3590  ENQRM=FNGRM 

3600  60  TO  3690 

3610  REM  STATEMENT  LABEL  11 

3620  FOR  11=1  TO  N 

3625  PRFACT  < 1 1 >=1 . 0 

3630  FOR  JJ=1  TO  N 

3640  A( I I , JJ>=H< I I , JJ> 

3650  REM  STATEMENT  LABEL  12 

3660  f'EXT  JJ 

3670  NEXT  II 

3680  EH0RM=1 . 0 

3690  REM  STATEMENT  LABEL  13 

3700  REM 

3710  SUBEND 


4006  SUB  HESQR<N, NM,  EUR<  > , EUI < > , SUBDI A<  > , INDICC  > , EPS, EX) 

4010  REM 
4020  REM 

4030  IF  N-2<0  THEN  4780 
4040  IF  N-2=0  THEN  4060 
4050  IF  N-2>0  THEN  4080 
4060  REM  STATEMENT  LABEL  1 
4070  SUBDIA<1>=A<2, 1> 

4080  60  TO  4730  ::  REM  60  TO  14 

4090  REM  STATEMENT  LABEL  2 

4100  M=N-2 

4110  FOR  KK=1  TO  M 

4120  L=KK+1 

4130  S=0. 0 

4140  FOR  I I I=L  TO  N 

4150  H(  1 1 1 , KK)=A< III, KK> 

4160  REM  STATEMENT  LABEL  3 
4170  S=S+ABS<A( III, KK) ) 

4175  NEXT  III 

4180  IF  SOABS<A<KK+l,KK)>THEN  4220  ::  REM  60  TO  4 
4190  SUBDIA<kK)=A<KK+l.KK> 

4200  H<KK+1,KK>=0.0 

4210  60  TO  4690  ::  REM  60  TO  12 

4220  REM  STATEMENT  LABEL  4 

4230  SR2=0. 0 

4240  FOR  I I I=L  TO  N 

4250  SR=A< III, KK) 

4260  SR=SR/S 

4270  A< 1 1 1 , KK)=SR 

4280  REM  STATEMENT  LABEL  5 

4290  SR2=SR2+SR*SR 

4300  NEXT  III 

4310  SR=SQR(SR2) 

4320  IF  A<L,KKK0.0  THEN  4340  ::  REM  60  TO  6 
4330  SR=—SR 

4340  REM  STATEMENT  LABEL  6 
4350  SR2=SR2-SR*A <L,KK) 

4360  A<L,KK)=A<L,KK)-SR 
4370  H<L,KK)=H<L,KK)-SR*S 
4380  SUBD I A ( KK ) =SR*S 
4390  X=S*SQR(SR2) 

4400  FOR  III-L  TO  N 
4410  H<III,KK)=H<III,KK)/X 
4420  REM  STATEMENT  LABEL  7 
4430  SUBO I A  < 1 1 1 ) =A  <1 1 1 , KK ) /SR2 
4440  NEXT  III 
4450  REM 

4460  REM  PREMULTIPLICATION  BV  THE  MATRIX  PR. 

4470  REM 

4480  FOR  JJJ=L  TO  N 
4490  SR=0. 0 
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4500  FOR  1 1 I=L  TO  N 

4510  SR=SR+A<  III, KK)*A< III, JJJ> 

4520  NEXT  III 
4530  FOR  1 1  I=L  TO  N 

4540  A<III,  JJJ  >=A<  III,  JJJ)-SUBDIA<  1 1 1  )*SR 
4550  NEXT  III 
4560  NEXT  JJJ 
4570  REM 

4580  REM  POSTMULTIPLICATION  BV  THE  MATRIX  PR. 

4590  REM 

4600  FOR  JJJ=1  TO  N 

4610  SR=0.0 

4620  FOR  I I I=L  TO  N 

4630  SR=SR+A< JJJ,  1 1 1  >*A<  III.JOO 

4640  NEXT  III 

4650  FOR  I II=L  TO  N 

4660  A< JJJ, 1 1 1 >=A< JJJ, 1 1 1 )-SUBDIA< 1 1 1 >*SR 
4670  NEXT  III 
4680  NEXT  JJJ 

4690  REM  STATEMENT  LABEL  12 

4700  NEXT  KK 

4710  FOR  KK=1  TO  M 

4720  A<KK+1,KK>=SUBDIA<KK> 

4730  NEXT  KK 

4740  REM  TRANSFER  OF  THE  UPPER  HALF  OF  THE  MATRIX  A  INTO  THE  ARRAY  H 

4759  REM  AND  THE  CALCULATION  OF  THE  SMALL  POSITIUE  NUMBER  EPS. 

4760  REM 

4770  SUBDIA<N-1 )=A<N, N-l ) 

4780  REM  STATEMENT  LABEL  14 
4790  EPS=0.0 
4800  FOR  KK=1  TO  N 
4810  INDIC(KK)=0 

4828  IF  KKON  THEN  EPS=EPS+SUBDIA<KKX2 
4830  FOR  III =KK  TO  N 
4840  H<KK, III >=A(KK, I I I > 

4350  REM  STATEMENT  LABEL  15 
4855  EPS=EPS+A  <  KK , 1 1 1 ) A2 
4860  NEXT  III 
4870  NEXT  KK 
4888  EPS=EX*SQR  <  EPS ) 

4890  REM 

4900  REM  THE  QR  ITERATIUE  PROCESS 
4910  REM 

4920  SHIFT=A<N,N-1> 

4930  IF  N<=2  THEN  SHIFT=0.0 

4940  IF  A<N,NX>0.0  THEN  SHIFT=0.0 

4950  IF  A<N-1,NK>0.0  THEN  SHIFT=0. 0 

4960  IF  A<N-1, N-l  >00.0  THEN  SHIFT=8.0 

4970  M=N 

4980  NS=0 

4990  MAXST=N*10 

5000  REM 
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5010  REM  TESTING  IF  THE  UPPER  HALF  OF  THE  MATRIX  IS  EQUAL  TO  ZERO. 

5820  REM  IF  IT  IS  EQUAL  TO  ZERO  THE  QR  PROCESS  IS  HOT  NECESSARY. 

5030  REM 

5040  FOR  I I 1=2  TO  N 
5050  FOR  KK=I I I  TO  N 

5060  IF  fi<III-l,KKK>0.0  THEN  5180  ::  REM  80  TO  18 

5070  NEXT  KK 

5088  NEXT  III 

5090  FOR  I I 1=1  TO  N 

5100  IND1CUII)=1 

5110  EUR<III>=A<III, 1 1 1 > 

5120  EUI<III)=0.0 
5130  NEXT  III 

5140  GO  TO  6620  ::  REM  GO  TO  37 
5150  REM 

5160  REM  START  THE  MAIN  LOOP  OF  THE  QR  PROCESS 
5170  REM 

5180  REM  STATEMENT  LABEL  18 
5190  K=M-1 
5200  M1=K 
5210  I=K 

5220  REM  FIND  THE  DECOMPOSITION  OF  THE  MATRIX 

5230  IF  K<0  THEN  6620  ::  REM  STATEMENT  37 

5248  IF  K=0  THEN  6309  ::  REM  STATEMENT  34 

5250  IF  K>0  THEN  5260  : :  REM  STATEMENT  19 

5260  REM  STATEMENT  LABEL  19 

5270  IF  ABS<A<M,K>X=EPS  THEN  6300  s:  REM  GO  TO  STATEMENT  34 
5280  IF  M-2=0.  THEN  6390  REM  GO  TO  35 
5290  REM  STATEMENT  LABEL  20 
5300  1=1-1 

5310  IF  ABS<A<K,I)K=EPS  THEN  5340  ::  REM  GO  TO  21 
5320  K=I 

5330  IF  K>1  THEN  5290  ::  REM  GO  TO  20 

5340  REM  STATEMENT  LABEL  21 

5350  IF  K=M1  THEN  6390  : :  REM  GO  TO  35 

5360  REM  TRANSFORMATION  OF  A  MATRIX  OF  THE  ORDER  GREATER  THAN  2 
5370  REM 

5380  S=A<M, M)+A<M1 , Ml J+SHIFT 

5390  SR=A<M, M)*A<M1 , Ml >-A<M, Ml >*A<M1 , M>+0. 25*SHIFT^2 
5400  A<K+2, K)-0. 0 
5410  REM 

5420  REM  CALCULATE  X1,V1,Z1  FOR  THE  SUBMATRIX  OBTAINED  BY  THE  DECOMPOSITION 
5430  REM 

5440  X=A<K, K>*<A<K,K)-S)+A<K,K+1 >*A<K+i , K>+SR 
5450  V=A<K+1,K>*<A<K,K>+A<K+1,K+1>-S> 

5460  R=ABS  <  X ) + ABS  <  Y ) 

5470  IF  R=0.8  THEN  SHIFT=A<M,M-1> 

5488  IF  R=8.0  THEN  5340  ::  REM  60  TO  21 

5498  Z=A<k>2,K+l>*A<K+l,)0 

5500  SHIFT=0. 0 

5510  NS=NS+ 1 

5520  REM 
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5530  REM  THE  LOOP  FOR  ONE  STEP  OF  THE  QR  PROCESS. 
5548  REM 

5550  FOR  1=K  TO  Ml 

5560  IF  I=K  THEN  5630  s:  REM  GO  TO  22 
5570  REM  CALCULATE  XR,YR,2R. 

5580  X=A< I, 1-1 > 

5598  V»A<I+1,I-1> 

5680  2=0.0 

5610  IF  I+2>M  THEN  5630  ::  REM  GO  TO  22 
5620  Z=A<I+2, I-l> 

5630  REM  STATEMENT  22 

5649  SR2=A8S  <  X  ) + A8S  <  V ) +ABS  <  Z  ) 

5650  IF  SR2=0.0  THEN  5690  ::  REM  GO  TO  23 
5660  X=X/SR2 

5678  V=V/SR2 
5680  Z=Z/SR2 

5698  REM  STATEMENT  LABEL  23 
5700  S=SQR<X*X+Y*Y+Z*Z> 

5710  IF  X<0.0  THEN  5730  ::  REM  60  TO  24 
5720  S=-S 

5730  REM  STATEMENT  LABEL  24 

5740  IF  I=K  THEN  5760  ::  REM  GO  TO  25 

5750  A<I, I-1)=S*SR2 

5760  REM  STATEMENT  LABEL  25 

5770  IF  SR2O0.0  THEN  5800  ::  REM  60  TO  26 

5780  IF  I+3>M  THEN  6230  : :  REM  GO  TO  33 

5790  GO  TO  6190  ::  REM  GO  TO  32 

5800  REM  STATEMENT  LABEL  26 

5810  SR=1.0-X/S 

5820  S=X-S 

5830  X=V/S 

5840  Y=Z'S 

5850  REM  PREMULTIPLY  BY  THE  MATRIX  PR. 

5860  REM 

5870  FOR  J=I  TO  M 

5880  S=A<I,J>+A<I+1, J>*X 

5890  IF  I+2>M  THEN  5910  ::  REM  GO  TO  27 

5900  S=S+A<I+2,J>*Y 

5910  REM  STATEMENT  LABEL  27 

5920  S=S*SR 

5930  A<I,J)=ACI,J>-S 

5940  ACI+l,  J)=A<I+1, J)-S*X 

5950  IF  I+2>M  THEN  5970  ::  REM  GO  TO  28 

5960  A<I+2,J)=ACI+2, J)-S*Y 

5970  REM  STATEMENT  28 

5980  NEXT  J 

5990  REM  POSTMULT I PLY  BY  THE  MATRIX  PR. 

6000  REM 
6810  L=I+2 

6020  IF  KMl  THEN  6040  ::  REM  GO  TO  29 
6030  L=M 

6040  REM  STATEMENT  LABEL  29 
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REM  60  TO  38 


6050  FOR  J=K  TO  L 
6060  S=A<J,I>+A<J,I+1>*X 
6070  IF  I+2>M  THEN  6090  :: 

6080  S=S+A<J, I+2>*V 
6090  REM  STATEMENT  LABEL  30 
6100  S=S*SR 
6110  ACJ, I>=A<J,  I>-S 
6120  A<J,I+i:>*A<:J,I+l>-S*X 
6138  IF  I+2>M  THEN  6150  ::  REM  GO  TO  31 
6140  A<J,I+2>*A<J,I+2>-S*V 
6150  REM  STATEMENT  31 
6160  NEXT  J 

6170  IF  I+3>M  THEN  6230  ::  REM  GO  TO  33 
6180  S=-A<I+3, I+2>*V*SR 
6196  REM  STATEMENT  32 
6200  A<I+3, I>*S 

6219  A<I+3,I+1>=S*X 

6220  ACI+3, I+2>=S*V+A<I+3,I+2> 

6230  REM  STATEMENT  33 

6240  NEXT  I 
6250  REM 

6260  IF  NS>MAXST  THEN  6620  ::  REM  60  TO  37 
6270  GO  TO  5188  ::  REM  GO  TO  18 
6280  REM 

6290  REM  COMPUTE  THE  LAST  El GENUAL UE 
6300  REM  STATEMENT  34 
6310  EUR<M>=A<M,M> 

6320  EUI<M>=0.0 
6338  INDIC<MJ>*1 
6340  M=K 

6350  GO  TO  5180  ::  REM  60  TO  18 
6360  REM 

6370  REM  COMPUTE  THE  EIGENUALUES  OF  THE  LAST  2X2  MATRIX  OBTAINED  BV  THE  DECOMPOS 

ITION 

6380  REM 

6390  REM  STATEMENT  LABEL  35 
6400  R=0.5*<A<K,IO+A<M.M>> 

6410  S=0.5*<A<M,M>-A<K,IO> 

6420  S=S*S+A<K,M)*A<M,K> 

6430  INDIC(K)=1 
6440  INDIC<M>=1 

6450  IF  S<0.0  THEN  6530  ::  REM  60  TO  36 
6460  T=SQR<S) 

6478  EUR(IO=R-T 
6480  EUR<M)=R+T 
6490  EUKK)=0.0 
6500  EUI<K>=0.0 
6510  M=M-2 

6520  GO  TO  5180  ::  REM  60  TO  18 
6530  REM  STATEMENT  36 


7000  SUB  REBLUEtN, NM, M, IUEC,fl< , > , UECRC , > , EUR< ) , EUI < > , IUORK< ) , WORK< ) , INDIC< > , EPS, 

ry  'I 

7010  REM 
7020  REM 

7030  UECR<1,IUEC)=1.0 

7040  IF  M=1  THEN  8340  : :  REM  GO  TO  24 

7050  EUALUE=EUR(IUEC) 

7060  IF  IVEC«M  THEN  7160  REM  GO  TO  2 

7070  K-IVEC+1 

7080  R=0.0 

7090  FOR  II=K  TO  M 

7100  IF  EUALU£OEUR<  II) THEN  7130  : :  REM  GO  TO  1 

7110  IF  EUI  Cl  1)00. 0  THEN  7130  ss  REM  GO  TO  1 

7120  R=R+3.0 

7130  REM  STATEMENT  1 

7140  NEXT  II 

7150  EUALUE=EUALUE+R*EX 

7160  REM  STATEMENT  2 

7170  FOR  K=1  TO  M 

7186  A<K,K)=A<K,K)-EUALUE 

7190  NEXT  K 

7200  REM 

7210  REM  GAUSSIAN  ELIMINATION  OF  UPPER-HESSENBER6  MATRIX  A 

7220  REM 

7230  K=M-1 

7248  FOR  11=1  TO  K 

7250  L=I 1+1 

7260  IUORK< I I )=0 

7270  IF  A<  1 1+1,11)08.0  THEN  7318  :  REM  GO  TO  4 

7280  IF  A<I 1, 11)00.0  THEN  7460  ::  REM  GO  TO  8 

7290  A<II,II)=EPS 

7300  GO  TO  7460  s :  REM  GO  TO  8 

7310  REM  STATEMENT  4 

7320  IF  ABS<A<II,II))>=ABS<A<II+1,II))THEN  7400  ::  REM  GO  TO  6 
7330  IWQRKX 1 1 )=1 
7340  FOR  >11  TO  h  ‘ 

7350  R=A<II,J) 

7360  A<II,J)=AUI+1,J) 

7370  REM  STATEMENT  5 

7380  A< I 1+1 , J)=R 

7390  NEXT  J 

7400  REM  STATEMENT  6 

7410  R=-A<II+1,IIVA<II,II) 

7420  A< I 1+1 , I I )=R 
7430  FOR  J=L  TO  M 
7440  A< I 1+1 , J)=A< I 1+1 , J)+R*A< I I , J) 

7450  NEXT  J 
7460  REM  STATEMENT  8 
7465  NEXT  II 
7470  REM 

7480  IF  A<M,M)O0.0  THEN  7520  ::  REM  GO  TO  9 
7490  A<M,M)=EPS 
7500  REM 
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7510  REM  THE  UECTOR  <1,1,... >  IS  STORED  IN  THE  PLACE  OF  THE  RIGHT  HAND  SIDE  COLU 
MN  UECTOR 

7526  REM  STATEMENT  9 
7530  FOR  11*1  TO  N 

7540  IF  I I>M  THEN  7570  ::  REM  GO  TO  10 

7550  UORK<II>=1.0 

7560  GO  TO  7590  ::  REM  GO  TO  11 

7570  REM  STATEMENT  10 

7580  UORK(II)*0.0 

7598  NEXT  II 

7600  REM 

7618  REM  THE  INUERSE  ITERATION  IS  PERFORMED  ON  THE  MATRIX  TIL  THE  NORM  ON  RHS  6T 

ER  THAN  BOUND 

7620  BOUND=0.01/<EX*N> 

7630  NS=0 
7640  ITER=1 
7650  REM 

7660  REM  THE  BACKSUBSTITUTION 
7670  REM 

7680  REM  STATEMENT  12 
7690  R=0.0 
7700  FOR  11=1  TO  M 
7710  J=M-I 1+1 
7720  S=WORK< J) 

7730  IF  J=M  THEN  7800  ::  REM  GO  TO  14 

7740  L=J+1 

7750  FOR  K=L  TO  M 

7760  SR=WORK<K> 

7770  REM  STATEMENT  13 
7780  S*S-SR*A<J,IO 
7790  NEXT  K 
7800  REM  STATEMENT  14 
7810  WORK< J)=S,'A< J,  J) 

7820  T=ABS<WORK<J>> 

7830  IF  R>=T  THEN  7850  ::  REM  GO  TO  15 
7840  R=T 

7850  REM  STATEMENT  15 
7855  NEXT  II 
7860  REM 

7870  REM  REM  THE  COMPUTATION  OF  THE  RIGHT-HAND  SIDE  UECTOR  FOR  THE  ITERATION  STE 
P. 

7880  REM 

7890  FOR  11=1  TO  M 
7900  WORK< I I )=WORK< I I )/R 
7910  NEXT  II 
7920  REM 

7930  REM  THE  COMPUTATION  OF  RESIDUALS 

7940  REM 

7950  R1=0.0 

7960  FOR  11=1  TO  M 

7970  T=0. 0 

7980  FOR  J=II  TO  M 

7990  T=T+A< I I , J>*UORK< J) 

8000  NEXT  J 


8010  T=ABS<T> 

8020  IF  R1>=T  THEN  8840  ::  REM  60  TO  18 
8030  R1=T 

8840  REM  STATEMENT  18 
8950  NEXT  II 

8860  IF  ITER=1  THEN  8880  ::  REM  60  TO  19 

8970  IF  PREUISOR1  THEN  8340  REM  60  TO  24 

8080  REM  STATEMENT  19 

8890  FOR  11=1  TO  M 

8100  REM  STATEMENT  20 

8110  OECR< I I , IUEC)=UORK( I I > 

8120  NEXT  II 
8138  PREUIS=R1 

8140  IF  NS=1  THEN  8340  ::  REM  60  TO  2454 
8158  IF  ITERM8  THEN  8360  ::  REM  60  TO  25 
8160  ITER=ITER+1 

8178  IF  R<B0UND  THEN  8200  ::  REM  60  TO  21 
8188  NS=1 

8190  REM  REM  GAUSSIAN  ELIMINATION  OF  THE  RIGHT-HAND  SIDE  UECTOR 

8200  REM  STATEMENT  21 

8210  K=M-1 

8220  FOR  11=1  TO  K 

8230  R=WORK( I 1+1 ) 

8240  IF  IWORK<II>=0  THEN  8280  :s  REM  GO  TO  22 
0250  WORK<II+l>=WORK<II>+WORK<II+l>*Aai+l,II> 

8260  WORK<II>=R 

8270  GO  TO  8380  ::  REM  GO  TO  23 
8288  REM  STATEMENT  22 

8290  WORK< 1 1+1 >=WORK< 1 1+1 >+WORK< 1 1 >*A< II+l,  II) 

8300  REM  STATEMENT  23 
8310  NEXT  II 

8328  GO  TO  7688  ::  REM  GO  TO  12 
8330  REM 

8340  REM  STATEMENT  24 

8350  INDIC< IUEC>=2 

8360  REM  STATEMENT  25 

8370  IF  M=N  THEN  8428  ::  REM  60  TO  27 

8380  J=M+1 

8390  FOR  II=J  TO  N 

8480  UECR<II,IUEC>=0.0 

8410  NEXT  II 

8428  REM  STATEMENT  27 

8438  SUBEND 
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9000  SUB  COMPUE<N,NM,M,  IUEC, A<,  ),UECR<,  ),H<,  ),EURO,EUIO.  INDICO,  IWORKO,SUBDIA 
(  > ,  WORKK  ) ,  U0RK2<  ) ,  WORK<  ) ,  EPS,  EX) 

9010  REM 

9020  REM  THIS  SUBROUTINE  FINDS  THE  IMAGINARV  EIGENVALUES  OF  THE  MAT IX 
9030  REM 

9040  FKSI=EVR<IUEC> 

9050  ETA=EUKIUEC) 

9060  IF  IUEC=M  THEN  9180  : :  REM  GO  TO  2 

9070  K=IU£C+1 

9080  R=0.0 

9090  FOR  I I-K  TO  M 

9100  IF  FKS I <>EUR< II) THEN  9130  ::  REM  GO  TO  1 

9110  IF  ABS<ETA)OABS<EUI  <  1 1 )  )THEN  9130  :  i  REM  GO  TO  1 

9120  R=R+3.0 

9130  REM  STATEMENT  1 

9140  NEXT  II 

9150  R=R*EX 

9160  FKSI=FKSI+R 

9170  ETA=ETA+R 

9180  REM  STATEMENT  2 

9190  R=FKSI*FKSI+ETA*ETA 

9200  S=2.0*FKSI 

9210  L=M-1 

9220  FOR  11=1  TO  M 

9230  FOR  J=II  TO  M 

9240  D=0. 0 

9250  A< J, I I )=0. 0 

9260  FOR  K=I I  TO  J 

9270  REM  STATEMENT  3 

9280  D=D+H(II,K)*H<K, J) 

9290  NEXT  K 

9300  REM  STATEMENT  4 

9310  A<II, J)=D-S*H<II, J) 

9320  NEXT  J 

9330  REM  STATEMENT  5 

9340  Aai,II)=A<II,II)+R 

9350  NEXT  II 

9360  FOR  11=1  TO  L 

9370  R=SUBDIA<II) 

9388  A<II+1, II)=-S*R 

9390  11=11+1 

9400  FOR  J=1  TO  II 

9410  A<J,II)=A<J,II)+R*H<J, II+l) 

9420  NEXT  J 

9430  IF  11=1  THEN  9450  :s  REM  60  TO  7 
9440  A<II+1, II-1)=R*SUBDIA<II-1) 

9450  REM  STATEMENT  7 

9460  FOR  J=II  TO  M 

9470  REM  STATEMENT  8 

9480  A<II+1, J)=A<II+1,J)+R*H<II, J) 

9490  NEXT  J 

9500  REM  STATEMENT  9 

9510  NEXT  II 
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9528  REM  REM  GAUSSIAN  ELIMINATION  ROUTINE 

9530  REM 

9540  K=M-1 

9550  FOR  11=1  TO  K 

9560  11=11+1 

9570  12=11+2 

9580  IWORK(II)=0 

9590  IF  II=K  THEN  9610  ::  REM  GO  TO  10 

9600  IF  A<  1 1+2, 11)00.0  THEN  9660  ::  REM  GO  TO  11 

9610  REM  STATEMENT  10 

9620  IF  A<  1 1+1,11)00.0  THEN  9660  ::  REM  60  TO  11 

9630  IF  A<II,  11)00.0  THEN  9950  ::  REM  GO  TO  18 

9640  A( I I , I I )=EPS 

9658  GO  TO  9958  : j  REM  GO  TO  18 

9660  REM  STATEMENT  11 

9670  IF  II=K  THEN  9730  ::  REM  GO  TO  12 

9680  IF  ABS<A<II+1, II))>=ABS<A<II+2,  II))THEN  9730  : :  REM  GO  TO  12 

9690  IF  ABS<A(II, II))>=ABS<A<II+2, II))THEN  9870  ::  REM  GO  TO  16 

9700  L=I 1+2 

9710  IW0RKXII)=2 

9720  GO  TO  9770  ::  REM  GO  TO  13 

9730  REM  STATEMENT  12 

9740  IF  ABS<A<II, II))>=ABS<A<II+1,  II))THEN  9840  ::  REM  GO  TO  15 

9750  L— I 1+1 

9768  IUORKC I I )=1 

9770  REM  STATEMENT  13 

9780  FOR  J=II  TO  M 

9790  R=A<II, J) 

9800  A<II, J)=A<L, J) 

9810  REM  STATEMENT  14 

9820  A<L,J)=R 

9830  NEXT  J 

9840  REM  STATEMENT  15 

9850  IF  1 1  OK  THEN  9870  ::  REM  GO  TO  16 

9860  12=11 

9870  REM  STATEMENT  16 
9880  FOR  L= 1 1  TO  12 
9890  R=-A<L,II)/A<II,II) 

9900  A<L, II)=R 
9910  FOR  J=I1  TO  M 
9920  REM  STATEMENT  17 
9930  A<L, J)=A<L, J)+R*A<II, J) 

9940  NEXT  J 

9945  NEXT  L 

9950  REM  STATEMENT  18 

9968  NEXT  II 

9970  IF  A<M,M)O0.0  THEN  9990  ::  REM  GO  TO  19 

9980  A<M,M)=EPS 

9990  REM  STATEMENT  19 

10000  FOR  11=1  TO  N 

10010  IF  I I>M  THEN  10050  ::  REM  GO  TO  20 
10028  UECR< 1 1 , IUEC)=1 . 0 
10030  UECROI,  IUEC-1)=1.0 
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10946  GO  TO  10089  ::  REM  GO  0  21 
10050  REM  STATEMENT  20 

10069  U£CR< 1 1 , 1U£C>=0. 0 

10070  UECRCI I, IUEC-1 >=0.0 
10080  REM  STATEMENT  21 
10090  NEXT  II 

10100  REM  REM  INUERSE  ITERATION  IS  PERFORMED 
10110  REM 

10120  BOUND=0.01/<EX*N> 

10130  NS=0 

10140  ITER=1 

10150  FOR  11=1  TO  M 

10160  WORKC 1 1 )=H< 1 1 , 1 1 >-FKSI 

10170  NEXT  II 

10188  REM 

10190  REM  SEQUENCE  OF  COMPLEX  UECTORS 
18200  REM 

10210  REM  STATEMENT  23 

10220  FOR  11=1  TO  M 

10230  D=WORK< I I >*UECR< I I , IUEC) 

10240  IF  11=1  THEN  10268  ::  REM  GO  TO  24 
10258  D=D+SUBDI A< I 1-1 >*UECR< I 1-1 , IUEC> 

10260  REM  STATEMENT  24 
10270  L=I 1+1 

18280  IF  L>M  THEN  10320  :  s  REM  GO  TO  26 

10290  FOR  K=L  TO  M 

10300  D=D+H< 1 I ,K)*UECR<K, IUEC) 

10310  NEXT  K 

10320  REM  STATEMENT  26 

10330  UECR< 1 1 , IVEC-1 )=D-ETft*UECR< 1 1 , IUEC-1 > 

10340  REM  STATEMENT  27 
10350  NEXT  II 
10360  REM 

10370  REM  GAUSSIAN  ELIMINATION  OF  RIGHT  HAND  SIDE  UECTOR 

10380  REM 

10390  K=M-1 

10400  FOR  11=1  TO  K 

10410  L=II+IWORK<II> 

10420  R=UECR<L, IUEC-1) 

1 0430  UECR < L , I UEC- 1 > =UECR <11, IUEC-1) 

10440  UECR< I I , IUEC-1 >=R 

10450  UECR<II+1, IUEC-1 )=UECR< I 1+1, IUEC-1 )+A< I 1+1 , II)*R 

10460  IF  I I=K  THEN  10480  ::  REM  GO  TO  28 

10478  UECR< I 1+2, IUEC-1 )=UECR< I 1+2, IUEC-1 )+A< I 1+2, I I >*R 

10480  REM  STATEMENT  28 

10490  NEXT  II 

10500  REM 

10510  REM  THE  COMPUTATION  OF  THE  REAL  PART 
10520  REM 

10530  FOR  11=1  TO  M 
10540  J=M-I 1+1 
10550  D=UECR<J, IUEC-1) 

10560  IF  J=M  THEN  19620  :: 


REM  GO  TO  38 
1 6 1 


10570  L=J+1 
10580  FOR  K=L  TO  M 
10590  D1=A<J,K) 

1 0600  0=0-0 1 *UECR <K, IUEC-1 ) 

10610  NEXT  K 

10620  REM  STATEMENT  36 

10630  UECRCJ, IUEC-1 )=D/A<J,J) 

10640  REM  STATEMENT  31 
18650  NEXT  II 
10660  REM 

10670  REM  THE  COMPUTATION  OF  THE  IMA6INARV  PART 
10680  REM 

10690  FOR  11=1  TO  M 

10700  D=WORK< I I )*UECR< I I , IUEC-1 > 

10710  IF  11=1  THEN  10730  s:  REM  80  TO  32 
10720  D=D+SUBDIA< 1 1 - 1 >*U£CR  C 1 1-1 , 1 UEC- 1 ) 

10730  REM  STATEMENT  32 
10740  L=I 1+1 

10750  IF  L>M  THEN  10796  ::  REM  80  TO  34 

10760  FOR  K=L  TO  M 

10770  D=D+H<II,K)*UECR<K, IUEC-1) 

10780  NEXT  K 

10790  REM  STATEMENT  34 

10800  UECR<  1 1 ,  IUEC)=<UECR<  1 1 ,  IUEO-DVETA 

10810  REM  STATEMENT  35 

10820  NEXT  II 

10830  REM 

10840  L=1 

10850  S=0.0 

10860  FOR  11=1  TO  M 

10870  R=UECR<  1 1 ,  IUEC)'“2+UECR<  1 1 ,  IUEC-1) ^2 
10880  IF  R<=S  THEN  10910  ::  REM  GO  TO  36 
10890  S=R 
10900  L=I I 

10910  REM  STATEMENT  36 

10920  NEXT  II 

10930  U=UECR<L, IUEC-1) 

10940  U=UECR<L, I UEC) 

10950  FOR  11=1  TO  M 
10960  B=UECR( I I , I UEC) 

10970  R=UECR< I I , IUEC-1 ) 

1 0980  UECR  < 1 1 , 1 UEC ) = <  R*U+B*U ) /S 
10990  UECR <11, IUEC-1 )=<B*U-R*U)/S 
11000  NEXT  II 
11010  REM 

11020  REM  THE  COMPUTATION  OF  THE  RESIDUALS 

11030  REM 

11040  B=8.0 

11050  FOR  11=1  TO  M 

1 I860  R=WORK(I I )*UECR< I I , IUEC-1 )-ETA*UECR< I I , IUEC) 
1 1070  U=WORK< 1 1 )*UECR< 1 1 , IUEC)+ETA*UECR< 1 1 , IUEC-1 ) 
11080  IF  11=1  THEN  11110  ::  REM  GO  TO  38 
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1 1090  R=R+SUBDIA< I 1-1 >*UECR< I 1-1 , IUEC-1 > 

11100  U=U+SUBDIA<II-1)*UECR<II-1, IUEC> 

11110  REM  STATEMENT  38 
11128  L=I 1+1 

11130  IF  L>M  THEN  11180  ::  REM  GO  TO  40 
11140  FOR  J=L  TO  M 
11150  R=R+H<II,J>*UECR<J,  IOEC-O 
11160  U=U+H<1I, J>*UECR<J, lUEC; 

11170  NEXT  J 

11180  REM  STATEMENT  40 

11190  U=R*R+U*U 

11200  IF  B>=U  THEN  11220  ::  REM  60  TO  41 
11210  B=U 

11220  REM  STATEMENT  41 
11230  NEXT  II 

11240  IF  ITER=1  THEN  11260  ::  REM  GO  TO  42 
11250  IF  PREUISOB  THEN  11400  ::  REM  GO  TO  44 
11260  REM  STATEMENT  42 
11270  FOR  11=1  TO  N 
11280  W0RK1<II>=UECR<II, IUEC) 

11290  REM  STATEMENT  43 

11300  U0RK2<II)=UECR<II, IUEC-1) 

11310  NEXT  II 
11320  PREUIS=B 

11330  IF  NS=1  THEN  11450  ::  REM  GO  TO  46 
11340  IF  ITER>6  THEN  11488  ::  REM  GO  TO  47 
11350  ITER=ITER+1 

11360  IF  BGUND>SQR < S ) THEN  10210  ::  REM  GO  TO  23 
11370  NS=1 

11380  GO  TO  10218  ::  REM  GO  TO  23 
11390  REM 

11400  REM  STATEMENT  44 

11410  FOR  11=1  TO  N 

11420  UECRCII,  IUEC>=UIGRK1<II> 

11438  OECR< II, IUEC-1 )=W0RK2< II) 

11440  NEXT  II 
11450  REM  STATEMENT  46 
11460  INDIC<IUEC-1>=2 
11470  INDICaOEC>=2 
11480  REM  STATEMENT  47 
11490  SUBEND 
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12008  SUB  I NPUTT ( ROW , COL , WORK  < ,  ) ,  I ERROR > 

12010  REM  INPUT  IN  A  MATRIX 
12020  REM 
12030  PRINT 

12040  INPUT  "ENTER  #ROWS,#COLS  OF  MATRIX-": ROW , COL 
12050  REM 
12060  REM 
12070  REM 

12080  FOR  1=1  TO  ROW 
12090  PRINT  “ROW  #",I 
12100  FOR  J=1  TO  COL 
12118  PRINT  J,  "  " 

12120  INPUT  UORK<I,J> 

12130  NEXT  J 
12140  NEXT  I 
12150  SUBEND 
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10  REM  PROGRAM  TO  CALCULATE  F  MATRIX 
15  REM 
20  REM  BV  A 
25  REM 

100  REM  4TH  ORDER  RUNGA-KIJTTA  ROUTINE  FOR  UP  TO  AN  STH  ORDER  SYSTEM 
110  REM 

120  CALL  CLEAR 
130  REM 

140  OPEN  #1 : ”PIO" 

150  REM 

160  DIM  X<8, 8) ,  XUIGRK<8, 8) ,  XNEW<8, 8) ,  AMAT C8,  S) 

170  DIM  FI  <8, 8>  ,F2<8,8> ,  F3<!8, 8> ,  F4<8,8) ,  F(8,8) 

180  DIM  K1<8,8),K2<8,8>,K3<8,8> .K4<8,8) 

185  DIM  JJ<8, 3> , OUT<8, 8) , fl<8,8> , B(8, 8), FW0RK1(8,8>,FW0RK2(8, 8) 

190  REM 

200  REM  INITIALIZE  ALL  UARIABLES 

210  REM 

220  ANGLE=0.0 

230  FOR  1=1  TO  8 

240  FOR  J=1  TO  8 

250  X(I,J)=0.0 

270  NEXT  J 

280  NEXT  I 

290  REM 

300  CALL  CLEAR 
310  REM 

320  REM  ENTER  THE  INPUT  PARAMETERS 

330  REM 

340  INPUT  "ENTER  #  OF  INCREMENTS": I NCR 
350  DELTANGLE=2*P I  / 1 NCR 
360  REM 

370  INPUT  "ENTER  ORDER  OF  SYSTEM" : ORDER 
372  INPUT  "ENTER  IN  B:":BE 

374  INPUT  "ENTER  IN  MU:": MU 

376  INPUT  "ENTER  IN  GAMMA: ": GAMMA 

378  INPUT  "ENTER  IN  PSQUARE: " : PSQUARE 

379  PRINT  "INPUT  J  MATRIX" 

380  CALL  INPUTKORDER, ORDER,  JJC,  >,IERROR> 

331  CALL  CLEAR  ::  PRINT  "ENTER  F<0>  MATRIX" 

382  CALL  I NPUTT <  ORDER , ORDER ,  X  < ,  > ,  I ERROR  > 

383  CALL  CLEAR 

384  REM 

385  OPEN  #2: "DSK. DAT ATWO.FFMATR IX ".SEQUENTIAL, DISPLAY  , UPDATE, FIXED 

386  PRINT  #2: INCR+1 

387  PRINT  #2: BE  ::  PRINT  #2: MU  ::  PRINT  #2: GAMMA  ::  PRINT  #2: PSQUARE 

388  PRINT  #2: JJ<1, 1>: :  PRINT  #2: JJ< 1 , 2>: :  PRINT  #2:JJ<2, 1>::  PRINT  #2:JJ<2,2> 

389  PRINT  #2: ANGLE  ::  PRINT  #2:X<1,1>::  PRINT  #2:Xd.2>::  PRINT  #2:X<2, 1)::  PRIN 
T  #2:X(2, 2) 
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3961  REM  MASTER  LOOP 
408  REM 

410  FOR  N=1  T.O  I  NCR  STEP  1 
420  DISPLAY  ATC2,4>:  ,,N=“.N 
430  REM 

440  REM  FIRST  RUNGA-KUTTA  CONSTANT 
450  REM 

460  CALL  FCNNCANGLE.XC, >,JJC, >,F1C, > , BE, MU, GAMMA, PSQUARE> 
470  REM 

480  FOR  1=1  TO  ORDER 
490  FOR  J=1  TO  ORDER 
500  K1CI,J)=DELTANGLE*F1C1,J> 

530  REM 

560  XWORKC I , J>=XC I , J)+0.5*K1 C I ,  J  ) 

570  NEXT  J 
580  NEXT  I 
590  REM 

600  REM  SECOND  RUNGA-KUTTA  CONSTANT 
610  REM 

620  CALL  FCNNCANGLE+0. 5*DEL TANGLE, XWORKC , > , JJC . > , F2C , ) , BE, M' 
630  REM 

640  FOR  1=1  TO  ORDER 
650  FOR  J=1  TO  ORDER- 
660  K2  C I , J  > =DELTAN6LE*F2 Cl, J) 

690  REM 

720  XWORKC I , J)=X< I , J>+0. 5*K2C I, J) 

730  NEXT  J 
740  NEXT  I 
750  REM 

760  REM  THIRD  RUNGA-KUTTA  CONSTANT 
770  REM 

780  CALL  FCNN C  ANGLE+0 . 5*DELTANGLE , XWORK  C , > , J J  C , > , F3  < , > , BE , f 
790  REM 

800  FOR  1=1  TO  ORDER 
810  FOR  J=1  TO  ORDER 
820  K3CI,J>=DELTANGLE*F3CI, J> 

850  REM 

880  XWORKC I, J>=XCI, J)+K3CI, J> 

890  NEXT  J 
900  NEXT  I 
910  REM 

920  REM  FOURTH  RUNGA-KUTTA  CONSTANT 
930  REM 

940  CALL  FCNNCANGLE+DELTAN6LE, XWORKC , ), JJC, >,F4C, >,BE,Ml 
950  REM 

960  FOR  1=1  TO  ORDER 
970  FOR  J=1  TO  ORDER 
980  K4CI,J>=DELTANGLE*F4CI, J> 

990  NEXT  J 
1000  NEXT  I 
1010  REM 
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1020  REM  EVALUATION  OF  STATE  VARIABLES 
1030  REM 

1040  FOR  1=1  TO  ORDER 
1050  FOR  J=1  TO  ORDER 

1060  XNEWU. J>=X<I, J)+a/6)*<Kia, J>+2*K2<I.J)+2*K3U,J>+K4 
1070  NEXT  J 
1080  NEXT  1 
1090  REM 

1100  REM  UPDATE  STATE  VARIABLES 
1110  REM 

1115  REM  PRINT  #1:"  ANGLE= " ,  < ANGLE+DELT  ANGLE >  *  1 80/'  <  P I ) 

1116  IF  (N/10>=INT<N/10)THEN  1117  ELSE  1120 

1117  PRINT  #2; ANGLE+DELTANGLE 
1120  FOR  1=1  TO  ORDER 

1130  FOR  J=1  TO  ORDER 
1140  Xa,J>=XNEW<I,  J> 

1145  REM  PRINT  #1,  USING  “  F<#,  :  I ,  J,  XNEUK I ,  J> 

1150  NEXT  J 
1160  NEXT  I 

1162  IF  <N/10>=INT(N/10>THEN  1165  ELSE  1170 

1165  PRINT  #2:X<1, 1 > : s  PRINT  #2:X<1,2>::  PRINT  #2:X<2,l):s  PRINT 
1170  REM 

1180  ANGLE=ANGLE+DELTANGLE 
1190  REM 
1200  NEXT  N 
1210  REM 

1220  REM  PRINT  OUTPUT 

1230  REM 

1231  PRINT  #1:" - 

___  II 

1232  PRINT  #1 

1233  PRINT  #1, USING  "  BE=##.##M~'  MU=##.##^A,,:BE,MU 

1234  PRINT  #1 

1235  PRINT  #1,  USING  "  GAMMA=## •  ##/s,/SAA  PSQUARE=##.##AAAA":  GAMMA 

1236  PRINT  #1 

1237  PRINT  #1 

1240  FOR  1=1  TO  ORDER 
1250  FOR  J=1  TO  ORDER 

1260  PRINT  #1, USING  “  x <#,#)=###. I,  J,X<I,  J) 

1270  NEXT  J 
1280  NEXT  I 
1290  PRINT  #1 
1300  PRINT  #1 
1310  PRINT  #1 
1320  REM 

1325  REM  NEXT  GAMMA 

1326  REM 
1330  CLOSE  #1 
1335  CLOSE  #2 
1340  REM 
1350  END 
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1360  SUB  FCNN<T,XWORKX, ), JJ<, )„F<, ),BE,MU,6ftMMft,PSGUARE) 
1370  REM 

1380  CALL  MATR I Xfl  <  T , RMfiT  < ,  > ,  BE ,  MU , GftMMft,  PSQURRE  > 

1390  REM 

1400  CALL  MfiTRIXMULT<ftMflT<, >,XWORKC, ),2,2,2,FW0RK1<, )) 

1401  REM 

1402  CALL  MfiTRIXMULT<XWORK< , >, JJ(, >,2,2,2,FWQRK2<, >> 

1403  REM 

1404  CALL  MfiTRIXftDD<FWORKl < ,  ),FW0R(<2(<  >f 2,2„F<,  >,-!> 


1406  REM 
1410  REM 
1420  REM 
1430  SUBEND 
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2000  SUB  MATRIXA<PSI , MATA< , > , BE f MU. GAMMA, PSGUARE ) 

2001  REM 

2002  REM  ALGORITHM  TO  CALCULATE  STATE  TRANSITION  MATRIX 

2003  REM  FOR  A  PERIODIC  FUNCTION 

2004  REM  IN  PARTICULAR-HELICOPTER  BLADE,  NO  FEEDBACK 

2005  REM  UPDATED  16  AUG  1984 

2006  REM 

200?  DEF  ARCS I N  <  22 ) =ATN  C  ZZ/SQR  < 1 -ZZ*ZZ ) ) 

2008  PSIT-PSI 

2810  IF  PSIT>=0.0  AND  PSIKPI  THEN  2020  ELSE  2012 

2012  IF  PSIT>=PI  AND  PSIT«PI+ARCSIN<BE/MU)  )THEN  2038  ELSE  2014 

2014  IF  PS I T>=  < 2*P I -ARCS I N < BE/MU ) > AND  PSIT<<2*PI )-. 000001  THEN  2030  ELSE  2016 

2016  IF  PSITX2*PI>-. 000001  THEN  PSIT=0.0  ::  SO  TO  2020  ::  ELSE  GO  TO  2018 

2018  GO  TO  2040 

2020  REM  PRINT  #1:,,2020" 

2022  K=  U  /3  )  *BEA3*MU*C0S  <  PS I T  > + a /4 ) *BE  ~2*MU^2*S  I N  <  2*PS  I T  ) 

2823  REM 

2024  C*  d  /4  >  *BE/'4+  a  /3  >  *BEA3*MU*S  I N  <  PS I T  > 

2025  REM 

2026  GO  TO  2050 

2030  REM  PRINT  #1:M2030" 

2031  REM 

2032  K=  a  /3 > *BEA3*MU*C0S <PSIT)+0/4) *BE~2*MUA2*S I N < 2*PS I T > +MU~4* <-a/12)*SINC2*P 
SIT)+< 1/24)*SIN<4*PSIT) ) 

2033  REM 

2034  C=< l/4)*BEA4+a/3)*BE~3*MU*SIN<PSIT)+MUA4*< < l/16)-< l/12)*C0SC2*PSIT)+( 1/48) 
*COS<4*PSIT)) 

2035  REM 

2036  GO  TO  2050 

2040  REM  PRINT  #1:"2040“ 

2041  REM 

2042  K=-  <l/3)  *BE~3*MU*C0S  <  PS  I T  )  -  <  1  /4  >  *BE/'2*MU''2*S  I N  <  2*PS  I T  ) 

2043  REM 

2044  C=-  <  1  /4  >  *BE A4-  0/3)  *BE/'3*MU*S  I N  <  PS  I T  ) 

2045  REM 

2058  REM  CONSTRUCT  MATRIX  A 

2051  REM 

2052  MATA<1,1)=0 

2053  MATA<1,2)=1 

2054  MATA  <  2 , 1 ) — <  GAMMA/2  )  *  <  2*PSQUARE/QAMMA+K  > 

2055  MATA  <  2 , 2 ) =- < GAMMA*C/2 ) 

2056  REM 

205?  REM  RETURN  TO  MAIN  PROGRAM 
2060  SUBEND 
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4000  SUB  INPUTT (ROW, COL, WORK( . ) , IERROR> 
4010  REM  INPUT  IN  ft  MfiTRIX  DwTh 
4020  REM 

S  K  "ENTER  »RWS,.COLS  OF  MATRIX-" 

4050  IF  R0W>8  THEN  4040 

4060  IF  ROWC2  THEN  I ERROR* 1 

4880  FOR  1*1  TO  ROW 

4090  PRINT  “ROW  #",I 

4180  FOR  J®1  TO  COL 

4110  PRINT  J,"  " 

At TWPliT  UORK<I.J> 


4130  NEXT  J 

4160  SUBEND  ::  REM  FROM  THE  SUBROUTINE 
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7000  SUB  MATRIXMULT< A< , > , B< , ) . ROWA. COLA. COLB. OUTc  , >  > 

7001  REM 

7002  REM  MATRIX  MULTIPLICATION  ROUTINE 

7003  REM  INPUTS  MATRIX  "A".  MATRIX  “B" . DIMENSION  OF  AAND  COLUMN  DIMENSION  Oc 
B 

7004  FOR  1=1  TO  ROWA 

7005  FOR  J=1  TO  COLA 

7006  OUTCI. J>=0.0 

7007  NEXT  J 

7008  NEXT  I 

7009  REM 

7010  FOR  1=1  TO  ROWA 
7020  FOR  K=1  TO  COLB 
7030  FOR  J=1  TO  COLA 
7040  SUM=A<I,J>*8<J,K> 

7050  OUT < I , K>=OUT < I , JO+SUM 
7060  NEXT  J 

7070  NEXT  K 
7080  NEXT  I 
7090  REM 

7999  SUBEND  ::  REM  RETURN  TO  MAIN  PROGRAM 


9800  SUB  MATRIXADDCAC,  >,8<.  >,RGWA.COLA,OUT<,  >,SI6NN> 
9010  REM 

9020  REM  MATRIX  ADD I T I ON/SUBTRACT 1 ON  ROUTINE 
9030  REM 

9640  REM  SI6NN:  +1-ADDITION, -1 -SUBTRACT ION 
9050  REM 

9060  FOR  1=1  TO  ROWS 
9070  FOR  J=1  TO  COLA 
9080  IF  SIGNN=-1  THEN  9100 
9090  OUT < I . J)=A< I . J)+B< I , J> 

9095  GO  TO  9110 

9100  OUT<I,J>=A<I.J>-Ba,J> 

9110  NEXT  J 
9120  NEXT  I 
9130  REM 
9140  SUBEND 
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100  REM  FOURIER  COEFFICIENT  ANAL VS IS  ROUTINE 
110  REM  WRITTEN  IN  TI  EXTENDED  BASIC 
120  REM 

125  CALL  CLEAR 

130  REM  PROGRAM  READS  TAPE#2  INFORMATION  FROM  F-MATRIX  PROGRAM 
140  REM  AND  DERIUES  18  FOURIER  COEFFICIENTS  AND  PRINTS  OUTPUT 
158  REM 

160  DIM  V<2,2>,A<19,2,2>,B<19,2,2> 

170  REM 

171  FOR  K=1  TO  19 

172  FOR  1=1  TO  2 

173  FOR  J=1  TO  2 

174  A<K.I,J>=8.0 

175  B<K, I,J>=0.0 

176  NEXT  j 

177  NEXT  I 

178  NEXT  K 

180  OPEN  #1: "PIO" 

190  OPEN  #2 : " DSK . DATATWO . FFMATR I X “ , SEQUENT I AL ,  D I SPL A V  , UPDATE, FIXED 
200  REM 

210  INPUT  #2: I NCR 

220  PRINT  #1: “  INCREMENT  OF  RUN  EQUALS  INCR-1? ”  STEPS" 

230  REM 

240  INPUT  #2: BE, MU, GAMMA , PSQUARE 
250  PRINT  #1 

260  PRINT  #1 : "  FOR  THIS  RUN;" 

270  PRINT  #1, USING  "  B=##.##  MU=##.##  GAMMA=## . ##  PSQUARE=##. ##": B 

E, MU, GAMMA, PSQUARE 
280  REM 

290  INPUT  #2: LAMBDA 1 , DUMP, DUMP, LAMBDA2 
300  PRINT  #1 

310  PRINT  #1:"  POINCARE  EXPONENTS" 

320  PRINT  #1 

33©  PRINT  #1,  USING  "  OMEGA  1  =###.  #####^^  0ME6A2=### . 

LAMBDA 1.LAMBDA2 
340  REM 
350  PRINT  #1 

360  PRINT  #1 s "  SVSTEM  OF  EIGENUECTORS" 

370  PRINT  #1 

380  REM 

390  REM 

400  IMAGE 

402  N=0.0 

410  INPUT  #2: ANGLE 

420  INPUT  #2:V<1,1),V<1,2),V<2,1),V<2,2> 

430  PRINT  #1, USING  400: V<1, 1),V<1,2) 

440  PRINT  #1, USING  400:V<2, 1),V<2,2) 

442  REM 
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445  CALL  F0UR1ER<AN6LE,V<, >, A<, , >.B(, „ >,N> 

446  REM 
450  REM 
460  REM 

470  FOR  N=1  TO  36 
480  INPUT  #2: X 

490  INPUT  #2:V(1, I>,V<1.2),V<2. 1>,V<2,2> 

500  CALL  FOURIER<X,V<, ),A<f . >,8<, . >,N> 

510  NEXT  N 
520  REM 
530  REM 

540  FOR  K=1  TO  19 
550  FOR  1=1  TO  2 
560  FOR  J=1  TO  2 

570  ACK, I,J>=<2/368>*A<K,  I,  J>*<10^3> 

580  B<K, I, J)=<2/360>*B<K, I,J)*(10/3) 

590  NEXT  J 
600  NEXT  I 
610  NEXT  K 
620  REM 

630  REM  PRINT  OUT  RESULTS 
640  FOR  N=1  TO  19 
650  PRINT  #1 

668  PRINT  #1, USING  "  A##  COEFFICIENT" : N-l 

678  PRINT  #1 

686  PRINT  #1, USING  400: A<N, 1, 1>,A<N, 1,2> 

690  PRINT  #1, USING  400:A<N,2, 1>, fl<N,2.2> 

708  PRINT  #1 

710  PRINT  #1, USING  "  B##  COEFFICIENT": N-l 

720  PRINT  #1 

730  PRINT  #1, USING  400: B<N, 1, 1>,B<N, 1,2) 

740  PRINT  #1, USING  400:B<N,2, 1>,B<N,2,2> 

750  NEXT  N 
760  REM 
770  REM 
780  PRINT  #1 
790  PRINT  #1 
800  PRINT  #1 
810  CLOSE  #1 
820  CLOSE  #2 
830  REM 


880  SUB  FOURIER<X,V(,>,ft<,,>JB(:,  ,  >,N> 

890  REN 

900  REN  FOURIER  COEFFICIENT  DETERN I NATION 
910  REN 

912  IF  N=0.  THEN  CONSTANT® 1  ELSE  60  TO  914 

913  60  TO  920 

914  IF  N=36  THEN  CONSTANT® 1  ELSE  60  TO  916 

915  60  TO  920 

916  IF  INT(N/-2)=N/2  THEN  C0NSTANT=2  ELSE  C0NSTANT=4 
918  REN 

920  FOR  KK=1  TO  19 
930  FOR  11=1  TO  2 
940  FOR  JJ=1  TO  2 

945  DISPLAV  AT<5, 5> : N, CONSTANT.  V< 1 . 1 > 

950  REN 

960  A < KK , 1 1 , J J ) =A < KK , 1 1 , J J ) +CONSTANT*V ( 1 1 ,  J  J  > *COS < (KK-1 >*X> 
970  REN 

980  BCKK, I I , JJ)=BCKK, I I , JJ>+CONSTANT*V< I I , JJ>*SINt (KK-l >*X> 

990  REN 

1000  NEXT  JJ 

1010  NEXT  II 

1020  NEXT  KK 

1030  REN 

1040  SUBEND 
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100  REM  FOURIER  COEFFICIENT  ANALYSIS  ROUTINE 
110  REM  WRITTEN  IN  TI  EXTENDED  BASIC 
120  REM 

125  CALL  CLEAR 

130  REM  PROGRAM  READS  TAPE#2  INFORMATION  FROM  F -MATRIX  PROGRAM 
140  REM  AND  DERIUES  18  FOURIER  COEFFICIENTS  AND  PRINTS  OUTPUT 

145  REM  FOR  GAIN  FOURIER  SERIES  OF  3  CBD  MATRICES 

146  REM 
150  REM 

160  DIM  V<2, 2> , AC2, 2) , B<2,  2>, 6(2, 2> 

161  DIM  AM<2,2,\BM<2,2),I1<2,2>,  I2<2.2>,0UT<2,2> 

162  DIM  A61 < 19, 2, 2) , AG2<!19, 2, 2> , AG3c!  19,2, 2> , BG1 <19, 2, 2) , BG2< 19, 2, 2) , BG3< 19, 2, 2)  , 

163  DIM  BK2,  1 ) , B2<2, 1> , B3<2.  1 ) , Gl(4, 4>  . G2C4,  4>  . 63<4,  4) 

164  DIM  VINU<4,4>, AF< 19,2, 2>,BF< 19,2, 2> 

165  REM 

170  REM 

171  FOR  K=1  TO  19 

172  FOR  1=1  TO  2 

173  FOR  J=1  TO  2 

174  A61<K, I, J>=0.0  ::  AG2<K, I , J)=0. 0  ::  A63<K, I, J>=B.0 

175  BG1<K, I, J)=0.0  ::  BG2<K, I, J>=0.0  ::  BG3<K, I , J)=@. 8 

176  NEXT  J 

177  NEXT  I 

178  NEXT  K 

180  OPEN  #lj"PIO" 

1 90  OPEN  #2 s " DSK . DAT AT WO . FFMATRIX", SEQUENT I AL , D I SPLAY  , UPDATE , FIXED 
200  REM 

202  B1<1, 1)=1. 

203  Bl <2, 1 >=0. 

204  B2<  1 , 1 >=0. 

205  B2<2, 1 >=1 . 

206  63(1, 1>=1. 

207  B3<2, 1)=1. 

208  REM 

210  INPUT  #2: I NCR 

220  PRINT  #1:"  INCREMENT  OF  RUN  EQUALS  "jINCR-l;"  STEPS" 

230  REM 

240  INPUT  #2 : BE , MU , GAMMA , PSQUARE 
250  PRINT  #1 

260  PRINT  #ls"  FOR  THIS  RUN;" 

270  PRINT  #1, USING  "  B=##.##  MU=##.##  GAMMA=## . ##  PSQUARE*##. ##" ; B 

E. MU, GAMMA, PSQUARE 
280  REM 

290  I NPUT  #2 : LAMBDA 1 , DUMP , DUMP , LAMBDA2 
300  PRINT  #1 

310  PRINT  #1;"  POINCARE  EXPONENTS" 


OMEGA 1  =### .  #####'”' 


0ME6A2®###. 


320  PRINT  #1 
330  PRINT  #1, USING  " 

“ : LAMBDA 1 , LAMBDA2 
340  REM 
350  PRINT  #1 

360  PRINT  #1 : "  SVSTEM  OF  EIGENVECTORS" 

370  PRINT  #1 
380  REM 
390  REM 

400  IMAGE  ###.####^^''  ###.##*#AAAA 

401  N=0.0 

402  FOR  1=1  TO  2 

403  FOR  J=1  TO  2 

404  61U,J>=0.0 

405  G2(I,J>=0.0 

406  63( I , J >=0. 0 

407  NEXT  J 

408  NEXT  I 

409  REM 

410  INPUT  #2: ANGLE 

420  INPUT  #2: YCl. 1 > , Y( 1 , 2> , YC2, 1 > , V(2, 2> 

430  PRINT  #1, USING  400: V( 1 , 1 >,  V(  1 , 2> 

440  PRINT  #1, USING  406: Y<2, 1 > , V<2, 2> 

441  REM 

442  CALL  INVERSE <2,  V(,  >,  VINU< .  ')  > 

443  CALL  MATRIXMULT (VINLK , ) , 61 ( , > , 2, 2, 1 , 61 ( , ) > 

444  CALL  MATRIXMULT (VINV< , ) , B2( , >,2,2,1, 62( , ) > 

445  CALL  MATRIXMULT(VINU( ,  >, B3C,  ), 2, 2, 1 . 63( ,  )  > 

446  REM 

447  CALL  FOIJRI ERC ANGLE, 61  (,  >,AG1C, ,  >,BG1(, ,  ),N> 

448  CALL  FOURIER (ANGLE, G2( , ) , AG2( ,  ,  ) , BG2< . ,  > , N> 

449  CALL  FOURIER (ANGLE, G3( , > , A63( , , > , B63( , , > , N> 

450  REM 
460  REM 

470  FOR  N=1  TO  36 
480  INPUT  #2:X 

490  INPUT  #2:V(1, 1 > , Y( 1 , 2) , V(2, 1>,V(2,2> 

491  REM 

492  CALL  INVERSE (2, V( , >, VINU( , > > 

493  CALL  MATRIXMULT(VINV(, >,B1(, >,2,2, 1,61 <,» 

494  CALL  MATRIXMULT ( YINV( , >,B2(, >,2, 2, 1,G2(, >> 

495  CALL  MATRIXMULT (VINV( , > , B3( , ) , 2, 2, 1 , G3( , ) > 

496  REM 

497  CALL  FOURIER(X, G1 ( , > , A61 ( , , > , BG1 ( , , > , N> 

498  CALL  FOUR I ER ( X , G2 ( , !> , AG2 ( , , > , BG2 ( , , > ,  N > 

499  CALL  FOURIER(X,  G3( , ') ,  A63( , ,  > ,  BG3( , ,  > ,  N> 

500  REM 
510  NEXT  N 
520  REM 
530  REM 
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FOR  K=1  TO  19 
FOR  1=1  TO  2 
FOR  J=1  TO  2 

06 KK,I,J>=( 2  - 360 > *06 1<K, I,  J > * <  1 023 •' 
flG2<K',  i’  J)=<22360>*062<K\.  I ,  J>*<  1023> 

063  < K , I ’ J > = < 22360 >  *063 <  K , I , J >  ■ *  < 1 023 > 
eG 1 < K . I . J > = < 22360 > *66 1<  K , I , J >  * < 1 023 > 

662 <K,  I ,  J>=< 2  '360 > *662 <  K , I , J  >  *  <1 023 > 

663  <K, I , J>=< 22360 > *663  < K , I , J >  *  < 1 023  > 

NEXT  J 
NEXT  I 
NEXT  K 
REM 

REM  PRINT  OUT  RESULTS 
IM06E  ###.####^'''A  ###.####'S~'A 

FOR  N=1  TO  19 
PRINT  #1 

PRINT  #1, USING  "  A##  COEFFICIENT":N-l 

PRINT  #1 

PRINT  #1, USING  635: 061 <N,  1, 1),062<N,  1,  1>,063<.N,  1. 1> 
PRINT  #1 , USING  635: 061 <N, 2, 1 > , 06 2<N, 2. 1 > , 063<N,  2,  1 > 
PRINT  #1 

PRINT  #1, USING  "  6##  COEFFICIENT" :N-1 

i  PRINT  #1 

i  PRINT  #1, USING  635:661 <N, 1, 1>,BG2<N, 1, 1>,B63<N, 1,1) 
i  PRINT  #1 ,  USING  635:  BG1<N,  2,  1 > 662<N,  2, 1  > ,  663 CN,  2,  1  > 
i  NEXT  N 
I  REM 
i  REM 

)  PRINT  #1 
i  PRINT  #1 
J  PRINT  #1 
i  CLOSE  #1 
i  CLOSE  #2 
i  REM 
i  END 
)  REM 

J  REM - 

H  REM 
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880  SUB  F0UR1ER<X,6<,>,AF<,,>,BF<:,,  >,N> 

890  REM 

900  REM  FOURIER  COEFFICIENT  DETERMINATION 
910  REM  . 

912  IF  N=0.  THEN  CONSTANTS  ELSE  GO  TO  914 

913  GO  TO  920 

914  IF  N=36  THEN  CONSTANTS  ELSE  60  TO  916 

915  GO  TO  920 

916  IF  INT<N/2>=N/2  THEN  CONSTANTS  ELSE  C0NSTANT=4 
918  REM 

920  FOR  KK=1  TO  19 
930  FOR  11=1  TO  2 
940  FOR  JJ=1  TO  2 
945  DISPLAV  AT<5,5>:N 
950  REM 

960  AF  <  KK , 1 1 , J J  >  =AF  <  KK , 1 1 , J J  >  +CONST ANT*G  < 1 1 , J J  > *COS <<KK-1>*X) 
970  REM 

980  BFCKK, 1 1 , JJ>=BF<KK, 1 1 , JJ>+CONSTANT*G< 1 1 , JJ)*SIN<  CKK-1 >*X> 

990  REM 

1000  NEXT  JJ 

1010  NEXT  II 

1020  NEXT  KK 

1030  REM 

1840  SUBEND 


5800  SUB  INUERSE<N2,A<,),B<,>> 

5001  REM  GAUSS- JORDAN  MATRIX  INUERSION  ROUTINE 

5010  REM 

5011  REM  IDENTIFIERS 

5013  REM  A  •  A  COEFFICIENT  MATRIX 

5014  REM  B  B  WORK  MATRIX 

5015  REM  B1  BIG  BIGGEST  UALUE 

5017  REM  D3  DETERM  DETERMINANT 

5018  REM  El  ERMES  ERROR  FLAG 

5018  REM  HI  HOLD  WORK  UARIABLE 

5028  REM  12  INDEX  WORK  MATRIX- 

5021  REM  13  I ROW  ROW  INDEX 

5022  REM  14  I COL  COLUMN  INDEX 

5023  REM  15  INURS  PRINT-INUERSE  FLAG 

5024  REM  N2  NCOL  NUMBER  OF  COLUMNS 

5025  REM  PI  PIUOT  PIUOT  INDEX 

5029  REM  END  OF  IDENTIFIERS 
5080  REM 

5890  E1=0  ::  REM  BECOMES  1  FOR  SINGULAR  MATRIX 
5100  15*0  ::  REM  PRINT  INUERSE  MATRIX  IF  ZERO 
5110  N3*l  ::  REM  NUMBER  OF  CONSTANT  UECTORS 
5120  FOR  1=1  TO  N2 

5130  FOR  J»1  TO  N2 

5140  B<I,J)=A<I,J) 

5150  NEXT  J 
5178  I2<I,3>=0. 

5180  NEXT  I 
5190  D3=l 

5200  FOR  1=1  TO  N2 
5210  REM 

5220  REM  SEARCH  FOR  LARGEST  PIUOT  ELEMENT 
5230  REM 
5240  B1=0. 

5250  FOR  J=1  TO  N2 

5260  IF  I2<J,3>=1  THEN  5350 

5278  FOR  K=1  TO  N2 

5288  IF  I2<K,3»1  THEN  6120 

5290  IF  I2CK, 3>=1  THEN  5340 

5300  IF  B1>=ABS<B<J,K))THEN  5340 

5310  I3=J 

5320  I4=K 

5330  B1=ABS<B<J,K>> 

5340  NEXT  K 
5350  NEXT  J 

5360  I2<I4,3)=I2U4,3>+1 
5370  I2<I, 1)=I3 
5380  12<l,2>=14 

5390  REM  INTERCHANGE  ROWS  TO  PUT  PIUOT  ON  DIAGONAL 

5400  IF  13=14  THEN  5540 

5410  D3=-D3 

5420  FOR  L=1  TO  N2 

5430  H1=B<I3,L> 


iV' 


5440  e<I3.L>=6<I4,L> 

5450  B<I4,L>=H1 
5460  NEXT  L 

5530  REM  DIM IDE  PIUOT  ROW  BV  PIUOT  ELEMENT 
5540  P1=B<I4. 14> 

5550  D3=D3*P1 
5560  BCI4,I4)=1 
5570  FOR  L=1  TO  N2 
5580  B<I4,L>aB<I4,L)/‘Pl 
5590  NEXT  L 
5640  REM 

5650  REM  REDUCE  NONPIUOT  ROWS 
5660  FOR  Ll  =  t  TO  N2 
5670  IF  Ll=14  THEN  5770 
5680  T*B<Ll,I4) 

5698  B(L1,I4>=0 

5700  FOR  L=1  TO  N2 

5710  B<L1,L>=B<L1,L>-B<I4,L>*T 

5728  NEXT  L 

5770  NEXT  LI 

5780  NEXT  I 

5790  REM 

5800  REM  INTERCHANGE  COLUMNS 
5810  REM 

5820  FOR  1=1  TO  N2 
5830  L=N2-I+1 

5848  IF  I2<L,1)=I2<L,2)THEN  5920 
5850  I3=I2(L, 1> 

5860  I4=I2<L,2) 

5870  FOR  K=1  TO  N2 
5880  H1=6<K, 13) 

5890  B<K, I3>=8<K, I4> 

5900  B<K, I4>=H1 

5918  NEXT  K 

5920  NEXT  I 

5930  FOR  K=1  TO  N2 

5940  IF  I2<K,3)01  THEN  6120 

5950  NEXT  K 

5960  E1=0. 

6800  GO  TO  6110 
6016  PRINT 

6020  PRINT  "  MATRIX  IMMERSION 

6038  FOR  1=1  TO  N2 

6040  FOR  J=1  TO  N2 

6850  PRINT  USING  M##.###~^";B<I,  J> 

6060  NEXT  J 
6070  PRINT 
6080  NEXT  I 
6085  BREAK 
6090  PRINT 

6100  PRINT  "DETERMINANT  =":D3 

6110  SUBEND  ::  REM  IF  INMERSE  IS  PRINTED 
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7000  SUB  MATRIXMULT<AM<, >,BM(, >, RQWA , COLA, COLB, OUT <, >) 

7001  REM 

74g  f^UTS  hE"B". DIMENSION  OF  »*D  COLO*,  DIMENSION  OF 

B 

7004  FOR  1*1  TO  ROUA 

7005  FOR  J=1  TO  COLA 

7006  OUT<I,J>*0.0 

7007  NEXT  J 

7008  NEXT  I 

7009  REM 

7010  FOR  1*1  TO  ROWfl 
7020  FOR  K=1  TO  COLB 
7030  FOR  J=1  TO  COLA 
7040  SUM=AM<I,J>*BM<J,K> 

7050  OUT (.  I ,  K)=OUT C I ,  IO+SUM 
7060  NEXT  J 

7070  NEXT  K 
7080  NEXT  I 
7090  REM 

7999  SUBEND  ::  REM  RETURN  TO  MAIN  PROGRAM 
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10  REM  PROGRAM  TO  CALCULATE  F  MATRIX 
15  REM 

20  REM  BV  A  FOURTH  ORDER  RUN6A-KUTTA 
25  REM 

100  REM  4TH  ORDER  RUNGA-KUTTA  ROUTINE  FOR  UP  TO  AN  8TH  ORDER  SVSTEM 
110  REM 

120  CALL  CLEAR 
130  REM 

140  OPEN  #1: "PIO" 

150  REM 

160  DIM  X(2,2>,XW0RK(2,2>,XNEW(2,2), AMAT (2,2> 

1 70  DIM  FK2,2>,F2<2,2\F3<2,2),F4<2,2),F<2.2> 

180  DIM  K1(8,8>,K2(8,8).K3(8,8>.K4(8,8> 

185  DIM  JJC , 2>,0UT(2,2>, A( 2,2). B( 2,2), FWORK 1(2,2), FWORK2( 2,2) 

186  DIM  FEEDBACK(2,2),XPRIME(2,2),XFU0RK(2,2).APRIME<2,2) 

188  DIM  AADD<2,2),BADD(2,2).AM(2,2),BM(2,2), 11(2,2), 12(2,2) 

189  DIM  XFEEDBACK(2, 2) , K1FB(2, 2) ,  K2FB(2, 2) . K3FB<2, 2) , K4FBC2, 2) 

199  DIM  FEED1(2,2),FEED2(2.2),FEED3(2,2),FEED4(2.2) 

191  DIM  FINUERSE(2,2) 

280  REM  INITIALIZE  ALL  UARIABLES 

210  REM 

220  ANGLE*0.0 

230  FOR  I«1  TO  2 

240  FOR  J»1  TO  2 

250  X<I,J>=0.0 

255  XFEEDBACK<I, J>=0.0 

256  IF  I*J  THEN  XFEEDBACKd,  J)=1.08 
270  NEXT  J 

280  NEXT  I 
290  REM 

300  CALL  CLEAR 
310  REM 

320  REM  ENTER  THE  INPUT  PARAMETERS 
330  REM 

340  INPUT  "ENTER  #  OF  INCREMENTS": I NCR 
350  DELTANGLE=2*P I ^ I NCR 
360  REM 

370  INPUT  "ENTER  ORDER  OF  SVSTEM": ORDER 
372  INPUT  "ENTER  IN  B:":BE 

374  INPUT  "ENTER  IN  MU:": MU 

376  INPUT  "ENTER  IN  GAMMA :": GAMMA 

378  INPUT  "ENTER  IN  PSQUftRE: " : PSQUARE 

379  PRINT  "INPUT  J  MATRIX" 

388  CALL  INPUTT (ORDER, ORDER, JJ<. >, IERROR) 

381  CALL  CLEAR  ::  PRINT  "ENTER  F<0>  MATRIX" 

382  CALL  INPUTT  (ORDER,  OR'>ER,  X< ,  > ,  IERROR) 

383  CALL  CLEAR 


384  REM 

385  REM  OPEN  #2: "DSK. DATATWO. FFMATRIX" , SEQUENTIAL. D I SPL A V  . UPDATE. FIXED 

386  REM  PRINT  #2:INCR+1 

38?  REM  PRINT- #2jBE  ::  PRINT  #2sMU  ::  PRINT  #2;GAMMA  ::  PRINT  #2:PSQUARE 

388  REM  PRINT  #2: JJ<1, 1>: :  PRINT  #2:JJ<1,2>::  PRINT  #2:JJ<2, 1>::  PRINT  #2:JJ<2,2 

) 

389  REM  PRINT  #2: ANGLE  ::  PRINT  #2s X< i . 1 >: s  PRINT  #2:X<l,2>:s  PRINT  #2:XC2,1>:: 
PRINT  #2:X<2,2> 

396  REM  MASTER  LOOP 
480  REM 

410  FOR  N=1  TO  I NCR  STEP  1 
420  DISPLAV  AT<2,4)j"N=",N 
430  REM 

440  REM  FIRST  RUNGA-KUTTA  CONSTANT 
450  REM 

460  CALL  FCNN< ANGLE, X< , > , XFEEDBACK< , >, JJ<, ),F1<, ), FEED1 ( , > , BE, MU, GAMMA, PSQUARE) 
4?0  REM 

480  FOR  1=1  TO  ORDER 
490  FOR  J=1  TO  ORDER 
500  K1 < I , J  >=DELTANGLE*F1 <  I ,  J  ) 

510  K1FB< I , J)=DELTANGLE*FEED1 < I , J) 

530  REM 

560  XWORKCI, J)=X<I,J)+0.5*K1<I,J) 

565  XFWORK< I , J ) =XFEEDBACK < I , J)+0. 5*K1FB< I , J) 

5?0  NEXT  J 
580  NEXT  I 
590  REM 

600  REM  SECOND  RUNGA-KUTTA  CONSTANT 
610  REM 

620  CALL  FCNN< ANGLE+0. 5*DELT ANGLE, XWORK< , ) , XFWORK< ,  > , JJ< , > , F2< , > , FEED2< , ) , BE, MU, 
GAMMA , PSQUARE ) 

630  REM 

640  FOR  1=1  TO  ORDER 
650  FOR  J=1  TO  ORDER 
660  K2< I , J)=DELTANGLE*F2< I , J) 

670  K2FB< I , J)=DELTANGLE*FEED2( I , J) 

690  REM 

720  XWORKC I , J)=X< I, J)+0. 5*K2< I , J? 

7 25  XFWORK (I, J) =XFEEDBACK < I , J ) +0 . 5*K2FB < I , J > 

7 30  NEXT  J 
740  NEXT  I 
750  REM 

760  REM  THIRD  RUNGA-KUTTA  CONSTANT 
778  REM 

780  CALL  FCNN< ANGLE +0. 5*DELT ANGLE, XUORKC , > , XFWORK < , > , , ) , F3< , ) , FEED3< , > , BE,  MU, 
GAMMA, PSQUARE) 

790  REM 

800  FOR  1=1  TO  ORDER 
810  FOR  J=1  TO  ORDER 
820  K3< I , J)=DELTANGLE*F3( I , J) 

825  K3FB< I , J ) =DEL T ANGLE*FEED3 < I , J) 

850  REM 
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880  XWORK<I, J>*X<I, J)+K3<I, J) 

885  XFWORKCI , J>=XFEEDBACK< I , J)+K3FB< I , J) 

890  NEXT  J 
960  NEXT  I 
910  REM 

920  REM  FOURTH  RUN6A-KUTTA  CONSTANT 
930  REM 

940  CALL  FCNN < ANGLE+DELTANGLE . XWORK  C . ) , XFUORK < . ) , J J < , ) , F4 < , ) , FEED4 < , > , BE , MU, 6AMM 
A , PSQUARE  > 

950  REM 

960  FOR  1=1  TO  ORDER 
970  FOR  J=1  TO  ORDER 
980  K4<  I , J>=DELTANGLE*F4< I , J ) 

985  K4FB< I , J>=DELTANGLE*FEED4 <1,  J) 

990  NEXT  J 
1000  NEXT  I 
1010  REM 

1020  REM  EUALUATION  OF  STATE  UARIABLES 
1030  REM 

1040  FOR  1=1  TO  ORDER 
1050  FOR  J=1  TO  ORDER 

1060  XNEWa,J>=X<I,J>+<l/6>*<Kia,J>+2*K2a,J>+2*K3a,J>+K4<I>J>> 

1065  XFEEDBACKCI , J >  =XFEEDBACK < I. J)+< 1^6>*<K1FB< I , J>+2*K2FB< I , J>+2*K3FB< I , J>+K4FB 

1070  NEXT  J 
1080  NEXT  I 
1098  REM 

1100  REM  UPDATE  STATE  UARIABLES 
1110  REM 

1115  REM  PRINT  #1:"  ANGLE=\  < ANGLE+DELTANGLE >*180/<P I > 

1116  IF  <N/10)=INT<N^10)THEN  1117  ELSE  1120 

1117  REM  PRINT  #2: ANGLE+DELTANGLE 
1120  FOR  1=1  TO  ORDER 

1130  FOR  J=1  TO  ORDER 
1140  X(I, J)=XNEW<I, J) 

1145  REM  PRINT  #1, USING  "  F<#, #>=###. I,  J,XNEW<I,J> 

1150  NEXT  J 
1160  NEXT  I 

1162  IF  <N/10>=INT<N/10)THEN  1165  ELSE  1170 

1165  REM  PRINT  #2:X<1,1>::  PRINT  #2:X<l,2^:s  PRINT  #2:X<2,1>::  PRINT  #2:X<2,2> 
1170  REM 

1180  ANGLE=ANGLE+DELT ANGLE 
1198  REM 
1200  NEXT  N 
1210  REM 

1220  REM  PRINT  OUTPUT 

1230  REM 

1231  PRINT  #1: “ - 

II 

1232  PRINT  #1 

1233  PRINT  #1, USING  "  BE=## .  MU=##.##~^"sBE,MU 

1234  PRINT  #1 


1235  PRINT  #1, USING  "  GAMMA=##.##A - 

1236  PRINT  #1 

123?  PRINT  #1 : "FEEDBACK  CONTROL  STYSTEM ' 

1238  PRINT  #1 

1248  FOR  1=1  TO  ORDER 

1250  FOR  J=1  TO  ORDER 

1260  PRINT  #1, USING  "  X<#, #>=###.###* 

1270  NEXT  J 

1288  NEXT  I 

1298  PRINT  #1 

1300  PRINT  #1 

1310  PRINT  #1 

1320  REM 

1325  REM  NEXT  GAMMA 

1326  REM 
1330  CLOSE  #1 
1335  REM  CLOSE  #2 
1340  REM 

1350  END 


PSGUARE=##.  :  GAMMA,  PSQUARE 


X  <#,#)=### .  1,3,  XFEEDBACK  <  I ,  J  > 


1360  SUB  FCNNCT, XWORKC , > , XFWORKC , > , JJ< ,  > ,  F< .  > .  XPRIMEC , > ,  BE. MU. GAMMA, PSQUARE > 
1370  REM 

1380  CALL  MATRIXACT,  AMATO.  BE,  MU,  GAMMA,  PSQUARE> 

1390  REM 

1400  CALL  MATRIXMULT <AMAT < , ) , XWORKX , >,2.2. 2, FU0RK1 < ,  >  ) 

1401  REM 

1402  CALL  MATRIXMULT <XWORK<,  ), JJC,  ),2,2,2,FWQRK2(.,  >> 

1403  REM 

1404  CALL  MATRIXADD<FW0RKlO,FW0RK2O,2,2,FO,-l> 

I486  REM 

1410  REM  FEEDBACK  CONTROL  SVSTEM  DETERMINED 
1420  REM 

1430  REM  CALL  INUERSE<2,XW0RK<, >, F I NUERSE < , >> 

1431  DET  2'XWORK  <  1 , 1 )  *XUORK  ( 2 , 2  >  -X  WORK  (.  1 , 2 )  *XWORK  <  2 , 1  > 

1432  FINUERSE<1,1>=XW0RK<2,2VDET 

1433  F I  NUERSE  <  1 , 2  >  =-XWORK  <1.2)  /'DET 

1434  FINUERSE<2,1>=-XW0RK<2.1VDET 

1435  FINVERSEC2, 2>=XW0RK< 1 , 1 VDET 
1440  REM 

1450  K=0. 109737 
1468  B1=0.0 
1470  82=1.0 
1480  REM 

1490  FEEDBACK (1,1) =K*B 1 *F I NUERSE  <  2 , 1  > 

1500  FEEDBACK < 1 , 2>=K*81*FINUERSE<2, 2> 

1510  FEEDBACKC2, 1 >=K*B2*FINUERSE<2, 1 ) 

1520  FEEDBACKS,  2>=K*B2*FINUERSE<2, 2> 

1530  REM 

1540  CALL  MATRIXADDCAMATC, ),FEED6ACK<, >,2,2, APRIMEC, >, 1> 

1550  REM 

1560  CALL  MATR I XMULT< APRIMEC, ),XFWORK<, ),2,2,2,XPRIME<, )> 

1570  REM 
1580  SUBEND 


2000  SUB  MATR  I  XA  <  PS  1 ,  MATA  < ,  > ,  BE ,  MU ,  GAMMA ,  PSGUARE  > 

2001  REM 

2002  REM  ALGORITHM  TO  CALCULATE  STATE  TRANSITION  MATRIX 

2003  REM  FOR  A  PERIODIC  FUNCTION 

2004  REM  IN  PARTICULAR-HELICOPTER  BLADE,  NO  FEEDBACK 

2005  REM  UPDATED  16  AUG  1984 

2006  REM 

2007  DEF  ARCSIN<ZZ>=ATN<2Z/SGR<1-ZZ*ZZ>> 

2008  PSIT=PSI 

2010  IF  PSIT>=0.0  AND  PSIKPI  THEN  2020  ELSE  2012 

2012  IF  PSIT>=PI  AND  PSI  T«PI+ARCSIN<  BE/MU)  >  THEN  2030  ELSE  2014 

2014  IF  PS I T>=  <  2*P  I  -  ARCS  I N  <  BE/MU  >  >  AND  PSIT«2*PI )-.  000001  THEN  2030  ELSE  2016 

2016  IF  PSITX2*PI>-. 000001  THEN  PSIT=0.0  ::  GO  TO  2020  ::  ELSE  GO  TO  2018 

2018  GO  TO  2040 

2020  REM  PRINT  #1:"2020‘’ 

2022  K=< 1/3) *BEA3*MU*C0S <  PS I T  > + <1 /4  > *BEA2*MUA2*S I N < 2*PS I T > 

2023  REM 

2024  C=  <1 /4  >  *BEA4+ <1  /3  >  *BEA3*MU*S I N  <  PS  I T  > 

2025  REM 

2026  GO  TO  2050 

2030  REM  PRINT  #1:"2030H 

2031  REM 

2032  K=  <  1  /3  >  *BE/'3*MU*C0S  <  PS I T  > + < 1 /4  >  *BEA2*MU''2*S  I N  <  2*PS  I T  >  +MU/S4*  < -  <1  /  1 2  >  *S  I N  <  2*P 
SIT>+<1/24>*SIN<4*PSIT>> 

2033  REM 

2034  C=<l/4>*€E'*4+<l/3>*BE''3*MU*SIN<PSIT>+MU'M*<<l/16>-<l/12:>*C0S<2*PSIT>+a/48> 
*C0S<4*PSIT>> 

2035  REM 

2036  GO  TO  2050 

2040  REM  PRINT  #1:"2040” 

2041  REM 

2042  K=-< l/3>*6E/'3*MU*C0S<PSIT>-<  1  /4 ) *BE''2*MU"'2*S I N < 2*PS I T ) 

2043  REM 

2044  C=-  <  1 /4  ■>  *BE A4-  <1/3  >  *BE'"3*MU*S  I N  <  PS  I T  ) 

2045  REM 

2050  REM  CONSTRUCT  MATRIX  A 

2051  REM 

2052  MATA<1,1>=0 

2053  MATA<1,2)=1 

2054  MATA<2, l>=-< GAMMA/2 >*<2*PSQUARE/GAMMA+K> 

2055  MATA  <  2 , 2 ) =- < GAMMA*C/2 ) 

2056  REM 

2057  REM  RETURN  TO  MAIN  PROGRAM 
2060  SUBEND 


4000  SUB  INPUTT (ROW, COL, WORK< , ) , I ERROR) 
4010  REM  INPUT  IN  A  MATRIX  DATA 
4020  REM 
4030  PRINT 

4040  INPUT  "ENTER  #ROU)S,#COLS  OF  MATRIX- 

4050  IF  R0W>8  THEN  4040 

4060  IF  RGW<2  THEN  I ERROR® 1 

4080  FOR  1=1  TO  ROW 

4090  PRINT  "ROW  #",I 

4100  FOR  J=1  TO  COL 

4110  PRINT  J,"  " 

4120  INPUT  WORK<I,J> 

4130  NEXT  J 
4150  NEXT  I 

4160  SUBEND  ::  REM  FROM  THE  SUBROUTINE 


5000  SUB  I NUERSE  <  N2 . A ( , > , B ( , > ) 

5001  REM  GAUSS- JORDAN  MATRIX  INUERSION  ROUTINE 


5010  REM 

5011  REM 

5013  REM  A 

5014  REM  B 

5015  REM  61 
501?  REM  D3 
5013  REM  El 

5019  REM  HI 

5020  REM  12 

5021  REM  13 

5022  REM  14 

5023  REM  15 

5024  REM  N2 

5025  REM  PI 

5029  REM 

5030  REM 
5090  E1=0  : 
5100  15=0  : 
5110  N3=l  : 


IDENTIFIERS 

fl  COEFFICIENT  MATRIX 

B  WORK  MATRIX 

BIG  BIGGEST  UALUE 
DETERM  DETERMINANT 
ERMES  ERROR  FLAG 
HOLD  WORK  UARIABLE 
INDEX  WORK  MATRIX 
I ROW  ROW  INDEX 

I COL  COLUMN  INDEX 

INURS  PR I NT- I NUERSE  FLAG 
NCOL  NUMBER  OF  COLUMNS 
PIUOT  PIUOT  INDEX 
END  OF  IDENTIFIERS 

:  REM  BECOMES  1  FOR  SINGULAR  MATRIX 
:  REM  PRINT  INUERSE  MATRIX  IF  ZERO 
:  REM  NUMBER  OF  CONSTANT  UECTORS 


5120  FOR  1=1  TO  N2 
5130  FOR  J=1  TO  N2 


5140  B(I,  J)=A<I,  J> 
5150  NEXT  J 
51?0  I2<I,3)=0. 
5180  NEXT  I 


5190  D3=l 

5200  FOR  1=1  TO  N2 
5210  REM 

5220  REM  SEARCH  FOR  LARGEST  PIUOT  ELEMENT 
5230  REM 
5240  B1=0. 

5250  FOR  J=1  TO  N2 

5260  IF  I2<J,3)=1  THEN  5350 

5270  FOR  K=1  TO  N2 

5280  IF  I2(K,3»1  THEN  6120 

5290  IF  I2(K,3>=1  THEN  5340 

5300  IF  B1>=ABS<BCJ,K))THEN  5340 

5310  I3=J 

5320  I4=K 

5330  B1=ABS(B<J.K)> 

5340  NEXT  K 
5350  NEXT  J 

5360  I2<14,3)=12<14,3>+1 
5370  12(1,1 >=I3 
53*30  12^»  I  2>-I4 

5390  REM  ’ INTERCHANGE  ROWS  TO  PUT  PIUOT  ON  DIAGONAL 

5400  IF  I 3=14  THEN  5540 

5410  D3=-D3 

5420  FOR  L=1  TO  N2 

5430  H1=B(I3,L) 

5440  B< 13, L)=B( 14, L) 
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5450  B<I4,L>=H1 
5460  NEXT  L  * 

5530  REM  DIVIDE  PIVOT  ROW  BV  PIVOT  ELEMENT 
5540  P1*B( 14, 14) 

5550  D3=D3*P1 
5560  B< 14, I4>=1 
5570  FOR  L=1  TO  N2 
5580  B<I4,L>=B<I4,LVP1 
5590  NEXT  L 
5640  REM 

5650  REM  REDUCE  NONPIVOT  ROWS 
5660  FOR  Ll=l  TO  N2 
5670  IF  L1=I4  THEN  5770 
568t*  T=BCL1,I4) 

5690  B<L1,I4)=0 

5700  FOR  L=1  TO  N2 

5710  B(L1 ,L)-B<L1 ,L>-B< I4,L>*T 

5720  NEXT  L 

5770  NEXT  LI 

5780  NEXT  I 

5790  REM 

5800  REM  INTERCHANGE  COLUMNS 
5810  REM 

5820  FOR  1=1  TO  N2 
5830  L=N2-I+1 

5840  IF  I2<L,1>=I2<L,2>THEN  5920 
5850  I3=I2<L, 1 ) 

5860  I4*I2<L,2) 

5870  FOR  K=1  TO  N2 
5880  H1=B<K, 13) 

5890  8<K, I3)=B<K, 14) 

5900  B(K, I4)=H1 

5910  NEXT  K 

5920  NEXT  I 

5930  FOR  K=1  TO  N2 

5940  IF  I2(K,3K>1  THEN  6120 

5950  NEXT  K 

5960  E1=0. 

6000  80  TO  6110 
6010  PRINT 

6020  PRINT  "  MATRIX  INVERSION 

6030  FOR  1=1  TO  N2 

6040  FOR  J=1  TO  N2 

6050  PRINT  USING  "##.###A~'~''":B<I,  J> 

6060  NEXT  J 
6070  PRINT 
6880  NEXT  I 
6085  BREAK 
6090  PRINT 

6100  PRINT  "DETERMINANT  ="jD3 

6110  SUBEND  ::  REM  IF  INVERSE  IS  PRINTED 


w'  L  .  '  V  „-.r- 


7000  SUB  MATRIXMULT <AM< ,  > ,  BM< , ) ,  ROWS,  COLA,  COLB,  OUT < ,  >  > 

7001  REM 

7002  REM  MATRIX  MULTIPLICATION  ROUTINE 

7803  REM  INPUTS  MATRIX  "A",  MATRIX  “B", DIMENSION  OF  AAND  COLUMN  DIMENSION  OF 
B 

7004  FOR  1=1  TO  ROUA 

7005  FOR  J=1  TO  COLA 

7006  OUT<I, J)=0.0 

7007  NEXT  J 

7008  NEXT  I 

7009  REM 

7010  FOR  1=1  TO  ROUA 
7020  FOR  K=1  TO  COLB 
7030  FOR  J=1  TO  COLA 
7040  SUM=AM< I , J)*BM< J,  tO 
7050  OUT < I ,  JO=OUT < I , JO+SUM 
7060  NEXT  J 

7070  NEXT  K 
7080  NEXT  I 
7090  REM 

7999  SUBEND  : :  REM  RETURN  TO  MAIN  PROGRAM 


9006  SUB  MATRIXADD<AADD< , > , BADD< , > , ROWA, COLA, OUTC , > , SI6NN> 
9010  REM 

9020  REM  MATRIX  ADD I T I ON^SUBTRACT I ON  ROUTINE 
9030  REM 

9040  REM  SIGNN:  + 1-ADDITION, -1 -SUBTRACT I ON 
9050  REM 

9060  FOR  1=1  TO  ROWA 
9070  FOR  J=1  TO  COLA* 

9080  IF  SIGNN=-1  THEN  9100 
9090  OUT < I , J)*AADD< I , J>+BADD< I , J> 

9095  GO  TO  9110 

9100  OUTCI, J>=AADD< I , J>-BADD< I , J> 

9110  NEXT  J 
9120  NEXT  I 
9130  REM 
9140  SUBEND 
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